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SUMMARY 


This  document  represents  an  interim  report  on  research  performed  under  Contract  No. 
F49620-89-C-0084  from  the  Air  Force  Office  of  Scientific  Research  to  Texas  A&M  University. 
The  period  covered  by  this  report  is  from  July  1,  1989  through  July  31,  1990.  The  research 
results  obtained  are  from  the  first  year  of  a  three  year  effort. 

- ^  Significant  progress  is  reported  on  analytical  and  computational  methodology  applicable  to 

dynamics  and  control  of  flexible  multibody  structures  structures.  Especially  significant  are  the 
following: 

'  -fr)  We  have  developed  new  analytical  and  numerical  results  pertaining  to  imposing  con¬ 
straints  in  multi-body  dynamical  modeling  and  numerical  simulation.  We  have  developed 
an  extension  of  existing  penalty  methods  for  constrained  multibody  dynamics,  including 
some  significant  convergence  proofs. ;  These  theoretical  developments  provide  a  rigorous 
foundation  for  existing  penal ty  jnethods  and  lead  directly  to  new  methods.  Of  special 
significance  are  the  analytical  and  numerical  studies  reported  for  linear  substructuring. 

—  <ii)~We  have  developed  a  power  principle  which  permits  the  efficient  construction  of 
stabilizing  control  laws  for  systems  described  by  nonlinear  systems  of  coupled  ordinary 
and  partial  differential  equations.^  It  is  not  necessary  to  first  spatially  discretize  the  partial 
differential  equations,  and  therefore  this  approach  is  immune  to  the  many  problems 
'  associated  with  truncation  and  spillover,  vis-a-vis  (for  example)  the  validity  of  the  stability 
boundaries.  Preliminary,  but  very  promising  results  have  been  obtained  using  this  prin¬ 
ciple,  including  analytical,  numerical  and  laboratory  experiments. 


1  (iii)  We  have  initiated  a  study  of  symbol  manipulation  methods  to  derive  polynomial-type 
nonlinear  feedback  control  laws  for  dynamical  systems  with  polynomial  nonlinearities.  A 
general  MACSYMA  symbolic  computer  code  has  been  developed  and  studies  are  under 
way  on  several  test  problems.  This  work  is  in  an  early  stage  of  development  and  the 
results  obtained  to  date  are  mixed;  we  have  found  several  attractive  control  laws,  but  for 
some  systems  surprising  and  as  yet  not  understood  failures  to  converge  have  occurred. 


Since  the  current  report  documents  the  interim  results  from  this  study,  we  anticipate  that  the 
above  results  will  have  been  significantly  extended  during  the  next  year. 


The  Investigators  for  this  effort  were  as  follows:  Professor  J.  L.  Junkins  served  as  Project . 


/ 


Director  and  Principal  Investigator.  Professor  A.  J.  Kurdila  and  Dr.  Z.  H.  Rahman  served  as 
Co-Principal  Investigators.  Three  Graduate  Research  Assistants  (GRAs)  participated  in  this 
project:  N.  Hecht,  R.  Menon,  and  S.  Hsu.  All  three  of  these  GRAs  are  pursuing  their  Ph.  D. 
dissertation  research  and  will  likely  be  near  completion  by  the  time  this  effort  is  completed. 


□ 


Organization  of  this  report  is  a  follows.  We  have  presented  the  detailed  technical  results  as 
attachments  to  this  report.  We  have  written  the  body  of  this  report  as  a  guided  tour;  following  a 
brief  introduction  in  Section  1,  the  sub-sections  of  Section  2  summarize  several  contributions 


with  reference  to  the  attachments.  Section  3  provides  concluding  remarks  and  discusses  some  ,/ 
promising  avenues  for  extending  the  research  discussed  in  this  report.  Codes 
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1.0  Introduction 


Dynamics  and  control  of  structures  research  has  recently  entered  an  exciting  era  characterized  by 
closer  coordination  of  theoretical,  computational,  and  experimental  work.  This  trend  is  at  least 
partially  in  response  to  the  challenge  to  apply  this  emerging  methodology  to  actual  space  struc¬ 
tures  within  the  next  decade.  A  key  to  accelerating  research  progress,  we  believe,  is  that  funda¬ 
mental  basic  research  needs  to  be  blended  with  some  constrained  degree  of  applied  research 
within  the  same  research  effort.  This  allows  some  carefully  designed  "reality  testing"  to  be 
carried  out  early,  providing  feedback  to  the  researchers  inventing  the  methodology,  and  to 
provide  "what  it  all  means"  insight  for  other  researchers  and  the  structural  dynamics  and  control 
community  at  large.  It  is  in  this  spirit  that  we  have  undertaken  the  present  research  project,  and 
therefore  we  iterate  between  basic  theoretical/ analytical  studies  and  carrying  out  constrained, 
simplified  applications  to  provide  some  benchmarks  for  evaluating  the  results.  In  this  interim 
report,  we  include  fundamental  analytical  results  as  well  as  some  of  the  numerical/experimental 
evaluation  studies  we  have  done  to  date.  The  technical  details  are  contained  in  the  several 
attachments. 

2.0  Technical  Accomplishments 

With  reference  to  the  attachments,  we  overview  the  research  contributions.  In  Section  2.1,  we 
discuss  some  significant  new  results  on  a  generalized  penalty  method  for  constrained  multibody 
dynamics.  Using  this  approach,  we  show  that  substructure  dynamical  models  can  be  married  into 
a  global  model  without  starting  over  and/or  the  necessity  of  inverting  large  linear  systems  to 
eliminate  lagrange  multipliers  or  solve  for  a  reduced  set  of  coordinates.  This  approach  has  some 
significant  advantages  over  existing  approaches.  New  analytical  convergence  proofs  and  conver¬ 
gence  bounds  are  presented  in  the  attachments.  In  Section  2.2,  we  summarize  results  fr  in  a 
method  for  designing  robust,  globally  stable  control  laws  for  nonlinear  distributed  parameter 
systems  using  a  power  principle.  We  show  how  to  use  this  method,  in  the  attachments,  to  design 
globally  stable  control  laws  for  near-minimum-rime  maneuvers  of  flexible  spacecraft  In  Section 
2.3,  we  overview  some  recent  work  we  have  done  on  an  approach  to  use  computer  symbol 
manipulation  to  derive  optimal  feedback  control  laws  for  nonlinear  dynamical  systems. 
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2.1  Lyapunov  Stable  Penalty  Methods  for  Multi-Body  Dynamics 


While  there  are  many  reasons  why  the  simulation  of  multiple  flexible  bodies  proves  such  a 
formidable  task  (the  existence  of  multiple  time  scales  for  rigid  and  flexible  degrees  of  freedom, 
high  dimensionality,  a  nonlinearly  varying  generalized  mass  matrix...)  one  fact  remains  after 
twenty  years  of  research:  no  completely  satisfactory  method  has  been  developed  that  efficiently 
and  reliably  accounts  for  the  "differential-algebraic"  nature  of  the  governing  equations.  It  is  well 
known  that  the  inclusion  of  constraints  using  (Langrange)  multipliers  can  be  the  cause  of  severe 
numerical  difficulties.  Researchers  in  numerical  analysis  have  investigated  general  integration 
schemes  applicable  to  a  wide  class  of  differential  algebraic  equations  [Petzold  1,2, 3,4],  Other 
researchers  in  the  field  of  multibody  simulation  have  suggested  specialized  means  of  detecting 
and/or  avoiding  these  inherent  conditioning  problems,  but  at  a  significant  computational  cost 
These  methods  include  the  nullspace  methods  described  in  [12],  the  range  space  methods  in  [21], 
Gears  method  [7],  a  particular  implementation  of  Kane’s  equations  [24]  and  Baumgarte’s  meth¬ 
ods  [3].  Some  of  the  most  novel  approaches  to  appear  recently  are  the  formulations  derived  by 
Park  [20]  and  [Bayo  1,2,3].  These  methods  use  a  penalty  approximation  of  the  Lagrange  multi¬ 
plier  terms. 

2.1.1  Basic  Ideas 

Penalty  methods  have  been  shown  to  be  theoretically  sound  for  certain  classes  of  boundary  value 
problems,  such  as  the  finite  deformation  of  an  incompressible  elastic  body  [19].  However, 
"typical"  hypotheses  involved  in  the  proofs  associated  with  these  such  problems  are  that 

(i)  the  governing  equations  are  obtained  from  the  Gateaux  differential  of  a  coer¬ 
cive  functional  (i.e.,  in  some  sense  the  functional  satisfies  a  growth  condition). 

(ii)  the  Langrange  multiplier  terms  satisfy  a  form  of  the  Babuska-Brezzi 

condition. 

It  is  not  difficult  to  cite  examples  in  multibody  formulations  of  dynamics  in  which  the  first 
criterion  is  not  satisfied.  But  one  should  not  conclude  that  penalty  methods  will  not  work  for 
multibody  dynamics  formulations,  only  that  the  existing  rigorous  tools  for  proving  existence  and 
convergence  of  the  method  cannot  be  applied  directly. 

The  fact  that  penalty  methods  must  be  employed  with  care  in  application  to  on  linear  formula¬ 
tions  of  dynamics  can  be  illuminated  by  a  few  simple  linear  examples.  Consider  first  the 
dynamics  of  the  free-free  finite  element  model  of  a  beam  shown  in  figure  (1).  The  usual  means 
of  obtaining  the  response  of  a  clamped-clamped  beam  would  be  to  use  the  constraint  conditions 
and  eliminate  the  redundant  degrees  of  freedom,  i.e.,  use  a  minimal  set  of  generalized  coor¬ 
dinates.  Alternatively,  one  can  retain  redundant  coordinates  and  use  Lagrange  multipliers  to 
account  for  the  constraint  forces  required  to  impose  the  constraints.  One  means  of  obtaining  an 
approximate  solution  to  the  problem  using  a  penalty  formulation  is  to  simply  adopt  the  estimate 
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as  is  commonly  employed  in  the  simulation  of  constrained,  coercive  variational  problems. 

This  procedure  can  lead  to  serious  numerical  difficulties  in  the  transient  problem  as  shown  in 
figures  (2)  and  (3).  In  this  collection  of  problems,  a  22-DOF  free-free  beam  shown  in  figure  (1) 
is  clamped  at  the  ends  using  the  simple  approximation  above.  In  the  first  figure,  the  constraint 
violation  at  the  ends  of  the  beam  is  plotted  as  a  function  of  the  penalty  parameter  e.  Clearly,  as 
the  penalty  parameter  approaches  zero,  the  constraint  violation  converges  to  zero  as  required. 
However,  figure  (3)  shows  that  this  approximation  of  the  multipliers  introduces  spurious  fre¬ 
quencies  that  become  unbounded  as  the  penalty  term  approaches  zero.  The  introduction  of  these 
extremely  high,  numerically  introduced  frequencies  makes  integration  of  the  approximate 
governing  equations  by  conditionally  stable  integration  schemes  impossible.  Just  as  importantly, 
an  eigenvalue  analysis  of  the  method  could  reveal  that  the  spurious  frequencies  are  inter-mixed 
with  the  true  structural  frequencies.  In  some  cases  it  is  difficult  to  determine  which  frequencies 
are  actual  structural  frequencies,  and  which  are  numerical  artifacts.  While  [Park]  and 
[Bayol,2,3]  have  provided  considerable  empirical  evidence  that  the  penalty  methods  can  be 
convergent  and  stable,  there  has  been  little  analytical  work  to  investigate  these  essential  features 
of  the  formulations. 


Throughout  this  research,  the  class  of  multibody  systems  under  consideration  consists  of  those 
that  can  be  represented  by  Lagrange’s  equations 


d  .dr.  ar  dv  „  a<t >/ 
d-t{Wk}~Wk+Wk=Qk  +  Wk 


h 


where  T  and  V  are  the  kinetic  and  potential  energies  of  the  system,  respectively,  qk  k=l..U  are 
the  generalized  coordinates,  Qk  are  generalized  forces  and  X,,  l  =  1...D  are  the  Lagrange  multi¬ 
pliers.  It  is  assumed  that  the  constraint  Jacobian  matrix 

atf)/ 

Wk 


has  been  derived  from  D  holonomic  constraints  having  the  form 


0/  (<7*)  =  0 


The  dynamical  system  of  Eqs.  (1),  (2),  whose  solution  is  the  true  motion  of  the  system,  is  re¬ 
ferred  to  as  the  "original,  dynamical  system." 

The  approach  presented  here  approximates  the  dynamics  of  the  original  system;  we  introduced 
penalized  potential  and  kinetic  energies  defined  by 
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Tt=Tt  (*,*)  =  r  +  i^p* 

Vt=Vc(qk)  =  V  +  j®T  a<t> 

where  a  and  (B  are  D  X  D  symmetric,  positive  definite  matrices  and  O  is  a  D  x  1  column  vector 
of  constraint  functions.  The  subscript  e  is  chosen  such  that  me  penalty  matrices  are 
parameterized  by  the  variable  e  >  0. 


a  =  a  (e) 

P  =  P  (e) 


The  penalty  matrices  are  most  simply  selected  to  be  diagonal  with 


Oo 

a‘J  ~  £  9/ 


where 

Oo  >0 

|3o>0 

One  can  also  introduce  a  "Rayleigh  constraint  penalty  dissipation  function"  via  the  definition 

Ft  =  j  <t>r  p.  <J> 


where  p  is  a  D  X  D  symmetric,  positive  definite  matrix 

|X  =  fx  (e) 

with  entries  ordinarily,  but  not  required  to  be,  defined  in  the  same  manner  as  a  and  P 


The  penalized  potential  and  kinetic  energies,  and  Rayleigh  dissipation  function  are  now  used  to 
generalize  Arnold’s  single  degree  of  freedom  developments  to  obtain  a  system  of  equations 
suitable  for  the  numerical  simulation  of  multiple  degree  of  freedom  multibody  systems. 


d  ( dTt  l  BTe  dV£  dF £ 
di{Wk}'Wk+Wk+Wk=Qk 


In  the  attachments  [7. 1-7.4]  it  is  shown  that  these  equations  arc  kinematically  equivalent  to  the 
constraint  violation  feedback  form  shown  below. 
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Notice  the  qualitative  affect  of  including  the  energy  penalty  terms  and  the  Rayleigh  constraint 
dissipation  function  is  to  invoke  "feedback  control"  generalized  forces  proportional  to  the 
constraint  violations  and  the  first  two  derivatives  thereof.  It  should  be  noted  that  the  above 
formulation  encompasses  the  starting  point  of  the  work  by  [Park],  as  well  as  the  stiffness,  damp¬ 
ing  and  inertial  methods  of  [Bayo]. 


We  introduce  a  generalization  of  Arnold’s  convergence  theorem: 


suppose 

(0 

(«) 


(«0 


(/v) 

(v) 


T(qk,qk)  =  T2{qk,qk) 

V  =  V(qk) 

qe(0)  =  q(0) 

<k  (0)  =  q(0) 

4>(q( 0),  $(0))  =  <D(<7  (0) )  = 

G  min  (a),  Gmin  (P)<->oo 


<35  e  <->  0 


then 


O') 

(«) 

m 


Et  =  E( 0)  =  rE  +  Ve  =  7(0)  +  V  (0) 
I<i>£l!2  +n  o, 

qc  <r->  q 


eii2  2  — ■-  — 77T7T  ++  <25  £  <— >  0 

fTlltl  (Ot),  fJmin  (p)  ) 


G5 


<t>F  +->  0 


where  of)  denotes  the  singular  values  of  (  ). 


In  addition  to  the  most  important  convergence  properties  above,  there  also  exist  a  wide  class  of 
systems  for  which  the  penalized  governing  equations  are  guaranteed  to  be  stable  in  the  sense  of 
Lyapunov.  Sufficient  conditions  for  the  stability  of  the  penalized  equations  are  summarized  in 
the  following  theorem  which  has  been  extracted  from  our  results  in  attachments  [7. 1-7.4].: 


> 

* 

♦ 

1 

t 

A 


Suppose 

(0  T(qk,qk)  =  T2(qk,qk) 

(»)  V=V(ft)  is  positive  definite  in  the  generalized  coordinates 

(iii)  6(0)  =  0(0)  =  0 

Then  a  sufficient  condition  that  the  penalized  governing  equations  are  stable 
in  the  sense  of  Lyapunov  is  that 

-<bTlL<t>  =  'LQkqk  <0 
k 

A  sufficient  condition  that  the  penalized  governing  equations  are  asymptotically 
stable  is  that  equality  holds  above  only  for 

q\  =qi  =  -  =  qk=o 

<7i  =  qi-  •••  =  qk  =0 


The  physical  interpretation  of  this  theorem  is  that  the  effects  of  the  constraint  penalty  is  to 
stabilize  the  open  loop  (Qk  =  0)  system,  since  ji  >  0. 

2.1.3  Imposing  Constraints  in  Linear  Substructuring  via  The  Penalty  Method 

The  above  general  results  can  be  sharpened  considerably  for  linear  systems.  In  particular, 
an  additional  theoretical  result  exploits  the  analogy  between  the  approximate  penalty  equations 
and  symmetric,  linear  quadratic  regulator  feedback  control.  Suppose  |3  =  0.  In  this  case,  the 
governing  linear  equations  can  be  written 

Mq  +  Cq  +  Kq  =  -  [^]  {^6  +  aO} 

Mq  +  Cq  +  Kq  =  - [^|]  ]  q  + 

By  inspection  the  penalty  terms  can  be  identified  as  symmetric  feedback  control. 

From  Joshi  [13],  the  symmetric,  linear  feedback  for  the  linear  system  above  is  the  optimal 
feedback  associated  with  the  following  performance  index: 
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or 


J  =  \  {XT WX  +  2 XTSU  +  UtRU}  dt 

0 


j  = 


W  S' 
S  R 


~X~ 

U 


dt 


where  X  is  the  2N  vector 


X  = 


q 

9 


With  the  definition  of  weight  matrices  shown  in  the  attachments  (7. 1-7.4),  the  performance  index 
becomes 


J  (2<bT [t<t>  +  2qT C q)  dt 

In  other  words,  when  P  =  0  and  the  system  is  undamped,  the  penalty  method  maximizes  a 
positive  measure  of  the  rate  at  which  the  (desired-to-be-zero)  constraint  energy  is  dissipated.  We 
can  also  use  LQR  methods  to  pursue  an  alternative  path  leading  to  the  stability  conclusions  of 
section  2.4.2. 

Frequency  domain  synthesis  of  substructure  models  is  also  possible  with  the  penalty 
approach.  The  undamped  eigenvalue  problem  associated  with  the  linear  system  is  shown  in 
attachments  [7. 1-7.4]  to  be 

[Mif]  «[w]]  +^[w+[w]  p[w]]]v:i0 

If  one  chooses 
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a  =  diagonal  (—  •••  -jf) 
P  =  diagonal  (—  •••  ~g) 


D  r  aoi7  r a<x>] 

p~l^j  LwJ 

p2  =p 


where  the  rows  of  the  Jacobian  matrix  have  been  orthononormalized  so  that  P  is  an  orthogonal 
projection  onto  the  space  of  admissable  configurations,  the  perturbed  eigenproblem  becomes 

{e(K  +  XtM)  +  (a0  +^p0)P}^£  =0 


By  introducing  the  abstract  angle  between  a  subspace  and  a  vector  as 


cos 


UPVtU 

11%  II 


an  error  bound  on  the  convergence  of  the  eigenproblem  is  derived  (attachments  [7. 1-7.4])  to  be 


CLp  +  K$o  I  COS  (P,%) 
OminiK  +  KM) 


The  application  of  the  approximate  eigenvalue  analysis  described  in  the  last  section  has 
been  applied  in  attachments  [7. 1-7.4]  to  a  two  substructure  system  having  120  degrees  of  free¬ 
dom  shown  in  figure  (4).  Figures  (5),  (6)  have  been  extracted  from  attachments  [7. 1-7.4]  and 
summarize  the  accuracy  of  the  spectral  estimates  for  a  range  of  penalty  parameters  varying  over 
eight  orders  of  magnitude.  The  first  figure  shows  that  the  norm  of  the  rms  relative  error  in  all 
108  natural  frequency  approximations  remains  below  10'3  for  all  values  of  e.  Figure  (6)  shows 
that  the  approximate  eigenvalues  obtained  using  LSPM  had  very  small  relative  errors  from  the 
true  frequencies.  These  errors  vary  from  roughly  3  or  4  digits  to  nearly  10  digits  of  accuracy  for 
different  values  of  the  penalty  parameter.  Thus,  the  method  is  extremely  robust  with  respect  to 
the  selection  of  the  penalty  parameter. 

2.1.4  Imposing  Constraints  in  Nonlinear  Multibody  Dynamics  via  the  Penalty  Method 

Several  simulations  have  shown  that  the  LSPM  are  extremely  effective  for  nonlinear 
multibody  systems.  Two  significant  features  of  the  method  make  it  especially  attractive  as  a 
means  of  constraint  stabilization. 
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Kir.lIRF.  (6)  NORMALIZED  ERROR  IN  EACH  FREQUENCY 


The  first  nonlinear  example  has  been  selected  to  show  that  the  penalty  method  can  be  an 
attractive  candidate  for  the  simulation  of  systems  having  configuration  dependent  singularities. 
Such  a  typical  system  is  the  closed  loop  planar  mechanism  shown  in  figure  (7).  Figure  (8)  shows 
the  time  histories  for  0,  calculated  by  the  range  space  method  and  penalty  method,  respectively. 
In  this  example  the  nonzero  constraint  violation  damping  term  has  been  the  selected  to  by  y  = 
60000.  As  is  evident  from  the  diagram,  the  penalty  method  remains  stable  throughout  the 
simulation,  while  the  range  space  method  eventually  diverges  at  a  singular  configuration.  The 
reason  for  the  improved  stability  when  the  penalty  method  is  employed  is  clear  from  the  con¬ 
straint  violation  plot  (9).  At  each  singular  configuration,  the  energy  transferred  to  the  constraint 
degrees  of  freedom  is  rapidly  dissipated,  so  that  the  norm  of  constraint  violation  returns  to  an 
acceptable  level.  It  should  be  noted  that  this  constraint  violation  stabilization  is  achieved  at  a 
cost:  as  the  constraint  energy  is  dissipated  after  each  singular  configuration  the  total  system 
energy  very  slowly  decreases,  instead  of  remaining  at  the  theoretically  constant  value.  The 
author  is  currently  investigating  the  use  of  energy-dissipation-rate-matching  integration  schemes, 
such  as  those  in  [5].  It  is  believed  that  this  combination  will  prove  to  be  a  powerful  tool  for 
configuration  singular  problems. 

As  a  general  observation,  it  should  be  mentioned  that  it  is  not  surprising  that  the  penalty 
formulation  described  herein  has  "regularizing"  characteristics;  it  resembles  the  method  of 
Tikhonov  regularization  employed  in  singularity-robust  pseudo-inverse  problems  [25]. 

The  utility  of  the  penalty  method  can  be  further  illustrated  in  application  to  nonlinear 
multibody  dynamics  problems.  Numerical  simulations  for  a  number  of  nonlinear,  natural, 
conservative  systems  have  verified  the  convergence  and  stability  theorems  derived  earlier. 
However,  the  assumptions  defined  in  the  theorem  above  are  somewhat  limiting  in  that  the 
convergence  is  guaranteed  only  for  natural,  conservative  systems.  Of  particular  interest  to  the 
authors  are  those  multibody  simulations  in  which  one  seeks  to  control  the  deployment  of  the 
system.  In  this  section  the  qualitative  behavior  of  the  above  numerical  procedure  for  the  simula¬ 
tion  of  dissipative  systems  is  considered. 

The  example  problem  is  shown  in  figure  (10)  is  discussed  further  in  attachments  [7. 1-7.4], 
The  essential  features  of  the  analysis  are  as  follows: 

(i)  The  convergence  and  stability  theorems  cited  earlier  guarantee  that  a  Lyapunov 
function  can  be  constructed  for  the  penalized  equations. 

(ii)  The  global  attractor  of  the  constraint  violation  trajectories  in  the  phase  plane  can 
be  constructed  using  the  Invariance  Principle. 

Simply  put,  this  theorem  asserts  that  the  trajectory  of  a  dynamical  system,  with  a  characterizing 
Lyapunov  function  defined  on  an  open  set  G  of  phase  space,  must  have  an  escape  time,  or 
approach  the  largest  positive  invariant  subset  M  of 
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FIGURE  (9)  NORM  OF  CONSTRAINT  VIOLATION 


Mc{(qk,qk)e  G|L  =  0} 

For  the  system  considered,  the  global  attractor  is  shown  in  the  attachments  [7. 1-7.4]  to  be  the 
family  of  ellipses  given  by 

.imp., 

K/a 

where  tc  is  defined  in  the  attachments  [7.1-7.4].  Figure  (1 1)  verifies  that  the  above  equation  does 
indeed  characterize  the  limit  cycle  to  which  the  constraint  violation  converges.  In  this  case,  a  = 
p  =  10000. 

Stronger  conclusions  can  be  obtained  for  this  problem  by  noting  that 

1 

and 


together  constitute  a  particular  solution  when  (see  attachments  [7. 1-7.4]) 

Kj  +  K2  =K 

The  different  limit  cycles  generated  by  specific  ratios  of  a  and  p  are  shown  in  the  attach¬ 
ments.  These  limit  cycles  obviously  exist  for  all  choices  of  a  and  P,  although  the  major  and 
minor  axes  may  be  sufficiently  small  for  a  particular  simulation  if  the  penalty  parameters  can  be 
made  sufficiently  large  without  encountering  numerical  difficulties.  Evidently,  if  an  initial  (or 
numerically  induced)  energy  spills  into  the  constraint  motion,  then  it  theoretically  remains  in  the 
invariant  subset  of  the  system.  Geometrically,  the  path  of  the  invariant  motions  for  this  example 
is  a  limit  cycle  corresponding  to  an  elliptical  translation  (parallel  displacement  of  the  bar)  without 
rotation.  Note  that  only  the  rotation  of  the  bar  is  controllable  by  the  feedback  torque;  the 
’’constraint  degrees  of  freedom,"  which  lie  in  the  invariant  subset,  are  uncontrollable  via  the 
dissipative  feedback  torque. 
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FIGURE  (10)  SIMPLE  3  DOF  SYSTEM  BEFORE  CONSTRAINT 


o.or-t-ooo  5.0E-004  i. nr  003  i.sf-003 


FIGURE  (11)  NORM  OF  ¥  LIMIT  CYCLE,  a  =  100,  (3  =  10000 


2.2  A  Power  Principle  for  Formulating  Feedback  Control  Laws 
for  Nonlinear  Distributed  Parameter  Systems 

2.2.1  Basic  Ideas 

Attachment  1  presents  a  novel  approach  we  have  developed  for  designing  robust,  globally  stable 
control  laws  based  upon  a  generalized  energy-rate  (power)  principle.  The  approach  is  more  akin 
to  classical  mechanics  than  control  theory,  however  it  elegantly  extends  the  Lyapunov  concepts 
used  for  stability  analysis  to  establish  a  control  law  design  method.  The  method  will  be  seen  to 
be  rather  democratic  in  that  it  applies  to  partial  as  well  as  ordinary  differential  equation  models, 
and  to  nonlinear  as  well  as  linear  models.  It  will  also  be  seen  to  lead  to  highly  robust  control 
laws  which  guarantee  stability  for  large  families  of  modeling  assumptions,  rather  than  simple 
parametric  or  additive  perturbations. 

Suppose  a  structural  system  is  acted  upon  by  m  actuators  which  can  impart  to-be-determined 
control  forces  (or  moments)  {Uj,  u2,  ...  ,  um]  ,  the  associated  coordinates  of  the  structural  loca¬ 
tions  at  which  the  actuator  acts  are  the  subset  of  system  (linear  or  angular)  coordinates  [qt,  q2,  ... 
,  qm).  The  first  step  in  the  development  requires  that  the  analyst  introduce  N+l  substructures  for 
the  purpose  of  describing  the  mechanical  system’s  energy  distribution  by  substructures.  The 
system’s  error  state  (current  position  and  velocity  state  minus  the  target  state)  is  defined  in  terms 
of  a  weighted  error  energy  function  having  the  form: 

N 

U  =  1  a, Ei  +aN.jf(qlt  q2,~qm>  Ci,  c2 . )  (2.2.1) 

i=0 

where  Ei  =  T  +  V,  is  the  total  mechanical  energy  of  the  ith  substructure,  Ti  is  the  ith  substruc¬ 
ture’s  total  kinetic  energy,  Vi  is  the  ith  substructure’s  total  potential  energy,  and  (a0,  aIt ... ,  aN+l) 
are  the  substructure  energy  weights;  these  weights  will  be  seen  below  to  parameterize  the  result¬ 
ing  feedback  control  laws.  The  function  f(qv  q2,  ...  ,  qm,  cv  c2,  ...)  is  introduced  for  generality, 
because  the  total  energy  is  not  always  positive-definite,  and  in  particular,  uncontrolled  bodies 
often  have  rigid  body  freedoms  which  must  be  eliminated  in  order  to  make  U  satisfy  the  funda¬ 
mental  necessary  condition  that  it  is  a  positive  definite  function  having  its  global  minimum  at  the 
target  motion.  Note  that  the  c,  denote  free  constants  which  may  be  chosen  subject  to  the  require¬ 
ment  that  /  remain  a  non-negative  function  ot  the  q’s.  In  the  present  form  of  Eq.  (2.2.1),  it  is 
implicitly  assumed  that  a  zero  energy  rest  state  is  the  target  motion,  however,  more  generally,  as 
developed  in  Attachment  1,  the  target  state  can  be  a  prescribed  reference  motion.  For  introduc¬ 
ing  the  ideas,  we  retain  the  simplest  form  of  Eq.  (2.2.1).  Consider  the  total  time  derivative  of  the 
error  energy  function,  this  is  obtained  by  differentiation  of  Eq.  (2.2.1)  to  obtain 


The  forces  (and  moments)  acting  on  the  ith  substructure  can  be  partitioned  into  four  subsets, 
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corresponding  to: 

(i)  Forces  having  a  potential  Vi ,  with  V  =  XV( , 

(ii)  Forces  which  do  no  work  on  the  ith  sub  structure  (e.  g.,  internal  forces), 

(iii)  Boundary  forces  and  moments  acting  on  the  ith  substructure,  and 

(iv)  The  control  force  acting  on  the  ith  substructure. 

Vectorially,  the  forces  acting  on  the  ith  substructure  can  be  written  as 

-H> 

F  i  —  V  V  i  +  F  non -working,  "l"  F  boundary  force  Si  BiU,  U  —COl{lli  Il2  ...  Mm}  (2.2.3) 

where  Bi  is  the  control  influence  matrix,  and  F boundary  force^  arises  from  interaction  with  adjacent 

substructures.  From  the  work-energy  principle  (applicable  with  the  same  generality  as  Newton’s 
Laws),  we  know  that  the  change  in  kinetic  energy  of  the  ith  substructure  is  given  by  the  work/ 
energy  equation  [26,  27,  28] 

7}  -  T0.  =worki=  Z-t  J  Fi  •  Ri  dt  (2.2.4) 

substructure-  * 
l  lo 


Substitution  of  Eq.  (2.2.3)  into  (2.2.4)  gives 

t 


E,  -  E0,  =  X 


substructure  • 

l 


t 


j  F boundary  force  • Ridt  +  J  BiU  *Ridt\  (2.2.5) 

| to  to 

where  £  =  T  +  Vt  is  the  total  energy  of  the  ith  substructure.  The  developments  of  this  section 
implicitly  assume  that  the  potential  energy  functions  are  non-negative,  for  the  more  general  case 
modifications  are  required.  Differentiation  of  Eq.  (2.2.5)  with  respect  to  time  gives  the  work-rate 
equation: 

dEi 


dt 


=  I  IF. 

substructure . 

I 


boundary  forces- 


i  •  Ri  +  Bi  U  • Rif  =  power 


(2.2.6) 


Note  that  the  above  discussion  can  be  easily  generalized  to  accommodate  boundary  and  control 
moments,  we  restrict  the  present  discussion  to  forces  for  simplicity.  A  more  general  discussion 
is  contained  in  the  attachments.  If  we  restrict  the  location  of  the  control  actuators  to  be  at  the 
boundaries  of  substructures,  then  for  a  natural  system  we  can  show  that  the  total  power  of  the  ith 
substructure  can  be  brought  to  the  form 
m 

^  =  Z  (2\y  +A‘iui)‘U 


j=  I 


(2.2.7) 


where  Qh  are  the  generalized  forces  associated  with  the  boundary  forces  and  moments,  and  A-u- 
is  the  jth  generalized  control  force  acting  on  the  ith  substructure.  The  qjs  are  displacement 
coordinates  of  the  actuator  locations.  Substitution  of  the  substructure  energy  rates  from  Eq. 
(2.2.7)  into  Eq.  (2.2.2)  provides  the  following  result  for  the  Liapunov  error  energy  rate  of 
change: 
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.  m  r  N  N  9f| 

0=%lSoa,Qbij  +£oa‘AiJ“J  +aN+iwJ\i,i 


(2.2.8) 


To  guarantee  that  Eq.  (2.2.8)  is  strictly  non-positive,  it  is  sufficient  to  require  that  each  of  the  m 
terms  in  the  [  ]  have  a  sign  opposite  to  ifj  .  The  simplest  choice  is  to  simply  set  each  of  the  [  ] 

terms  to  a  negative  constant  (-  kj)  times  qj  ,  this  gives  u  ~  kjqf  and  results  in  the  follow¬ 
ing  system  of  algebraic  equations 
N  N 

y  n  •  -4-  V/7-/1..I/.  -L  n . .  — ~  —  —  Jr  i  h  . 


SoaiQhij  +  /5)^'Ay UJ  +  Qn+1^j 


(2.2.9) 


These  m  algebraic  equations  can  be  solved  explicitly  for  the  m  control  functions  and  we  can 
establish  a  stabilizing,  constant  gain,  output  feedback  control  law  as 


Uj  =- 


(iy  + ^ + ka‘Q>,l)-i='’2 . m  (2-2,o) 


Note  that  the  above  control  law  feeds  back  measurements  of  position  (through  the  dependence  of 
/  on  the  q’s  ),  velocity,  and  boundary  forces  Qb  .  Notice  that  the  constant  gains  are 

V 

parameterized  as  a  function  of  the  energy  weights  (a,  )  and  the  constants  (ky  )  and  (cy  ).  Of 
course,  the  freedom  to  feed  back  boundary  forces  is  only  an  advantage  if  these  can  be  measured 
or  estimated  (via  strain  gauges  or  load  cells,  for  example).  Observe  that  the  boundary  forces  on 
adjacent  substructures  will  have  equal  magnitude  and  opposite  sign,  resulting  in  the  last  sum  in 


the  bracket  combining  in  pairs;  for  the  special  case  of  a  "chain  configuration",  for  example,  the 
N 

sum  has  the  structure:  £  a.Qb  _  (fl0-flj  )Qb  +  (flj  -  a2)  Qb  .+  ...>  ^us,  «  is  possible 
i — o  ty  .  j 


to  place  constraints  on  the  substructure’s  energy  weights  (in  this  case,  aimJ  -  a ()  to  either  allow  or 


eliminate  selected  boundary  force  feedback  terms;  note  that  the  details  of  this  summation  is  a 
function  of  the  system  topology  and  must  be  carried  out  specifically  for  each  application.  If  it 
can  be  shown  for  the  particular  system  that  the  generalized  actuator  location  velocities 

{<7l  >  <72  >•••»  }  vanish  identically  only  at  the  target  state,  then  it  is  evident  that  0  =  -  ^k,qf 

is  globally  negative  and  we  therefore  have  global  asymptotic  stability  of  the  closed  loop  system. 


The  above  discussion  may  appear  a  little  abstract,  but  as  is  evident  in  the  presentation  of  recent 
results  below,  it  can  be  applied  in  a  fairly  straightforward  way  to  achieve  some  very  attractive 
control  laws. 


2.2.2  Discussion  of  Recent  Results 

With  reference  to  figures  2.1  and  2.2,  we  consider  specializing  the  above  results  for  a  particular 
multi-body  maneuver  problem.  The  nine  body  configuration  and  modeling  assumptions  are 
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2.2  Solution  of  Problem  I 

We  seek  an  optimal  feedback  control  law  to  compute  u( x)  for 

x  =  - x  +  a x2  +  u  ,  a  >0 


( 2.2.29 ) 


if 

which  minimizes  the  performance  measure  ^  ~  J  j  C*2  +  W2  )  <i/ 

o 

The  Hamiltonian  for  this  system  is 

//=  i  (x2  +«2)  +  ^(-jc  +  ajt2+u) 


(2.2.30) 


(2.2.31) 


The  Pontryagin  necessary  conditions  are: 


M  =  -  X 


X  —  — g—  —  —x  +  X  -  2ax  X 


VII  .  2 

x  =  -yr  =  -x  -  X  +  a  xr 


(2.2.32) 


We  seek  a  feedback  form  of  the  control  law,  specifically,  we  seek  the  polynomial  gains  Ki  in 


-  u  =  X  =  I  Kpf 

i  =  1 


(2.2.33) 


Substituting  Eq.  (2.2.33)  and  its  time  derivative  into  Eqs.  (2.2.32),  we  find  the  homogeneous 
condition 

[tf,  -  2KX  -  K\  +  llx  +  [k2  -3(1+K,  )K2  +3<x/Ci  Jx2  + 

[^3  -4(1+*!  )Kj  +4ctK\  -  2K\)^  +  ...  +  [kN  -{N+l)(l+K\  )KN  -  FN(aJCr . KN.X  )]x»  =  0 

Since  the  above  equation  must  hold  at  every  point  in  the  state  space  (i.  e. ,  for  all  x),  we  conclude 
that  all  [  ]’ed  coefficients  must  vanish  independently,  this  provides  the  following  equations: 

ki  -2KX  -K\  +  1=0  =>  Ki 

k2  -3(1+*:,  )K2  =  -  3aK,  =>  K2 

k3  -4(l+Ki)Ki  =-4aKi  +2K\  =>  AT3  (2.2.34) 


kN  -(N+W+Ki  )Kn  =  Fn(  a,  Kx,K2,  ...JCN-i )  =>  K» 

Notice  the  following  structure  and  properties  of  Eqs.  (2.2.34)  satisfied  by  the  optimal  gains: 
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These  equations  can  be  solved  sequentially  for  the  Kt,  all  but  the  first  are  linear  equations. 
The  equations  can  solved  to  arbitrary  order,  since  we  have  developed  explicit  recursions 


for  Fn. 


®  The  first  equation  for  Kt,  is  a  scalar  Riccati  equation  (no  surprise  here!). 

®  If  we  impose  a  =  0  and  the  boundary  conditions  Kfc)  =  0  [consistent  with  X.(tf)  =0  as  a 

transversality  condition  for  x(t{)  free],  we  see  that  all  nonlinear  gains  vanish  identically  and 
therefore  the  above  is  indeed  a  direct  generalization  of  the  classical  linear  regulator  with  a 
quadratic  performance  index,  the  optimal  control  of  Eq.  (2.2.33)  reduces  to  u  =  -  Kt  x. 

®  If  tf  — >  ,  analogous  to  the  classical  steady  state  regulator,  we  can  show  that  all  Kt  ap¬ 

proach  constants  and  Eqs.  (2.2.34)  therefore  reduce  to  a  sequence  of  algebraic  equations 
for  the  constant  gains  (we  can  show  that  the  positive  real  root  is  the  proper  selection  for 
Kj),  and  the  solution  for  the  higher  gains  involves  only  simple  algebraic  operations.  For 
this  case,  with  a  =  1,  the  first  few  gains  are  numerically: 

KI  ~  0.41321,  K2  =  0.29289,  ^  =  0.17678,  K4  =  0.088388,  Ks  =  0.033145. 


Shown  below  in  Figure  2.3  is  the  performance  index  versus  the  order  N  of  the  feedback  control 
and  the  trajectories  of  x(t)  and  u(t)  for  typical  initial  conditions  [jc(0)=1.3].  As  is  qualitatively 
evident,  the  nonlinear  terms  are  constructive  and  convergence  is  rapid.  We  have  shown  that  the 
nonlinear  controls  for  N>3  are  globally  stable,  whereas  large  but  finite  domains  of  stability  are 
associated  with  the  linear  and  quadratic  feedback  control  laws.  This  solution  has  been  verified 
by  the  symbolic  manipulator  code.  As  we  indicate  below,  the  structure  of  this  simplest  scalar 
example  fully  generalizes,  so  that  the  above  observations  apply  in  a  much  more  general  context. 
However,  the  stability  and  convergence  issues  are  problem-dependent,  as  should  be  expected. 
We  have  established  a  means  to  generate  the  matrix  equivalent  of  the  sequence  of  scalar  equa¬ 
tions  ( Eqs.  (2.2.34  )  and  their  solutions  for  the  optimal  nonlinear  feedback  gains 

Solution  of  Problem  II 

We  seek  an  optimal  feedback  control  law  to  compute  u,(xIt  x2),  u2(xv  x2)  for  the  system: 
xi  =-xj  +X1X2  +  *2  +  My 

2  (2.2.35) 

JC2  =  -  X2  +  Xj  X2  +  X]  +  «2 

which  minimizes 

J  =  j  j  Z(x?  +uj)dt  =  \{  (fQx  +  urRu)dt  (2.2.36) 

0  i- 1  0 

Eqs.  (2.2.35)  can  be  written  as  i  =A\X\  +A2X2  +  Bu  ,  Q-R-l  (2.2.37) 
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Figure  2.2  Geometric/Kinematic  Notations  for  the  Equations  of  Motion 


S^stemLierangian: 

L  =  T  -V  -  Tkub+  T^-V^  =  +  -  '  / £ /, 


A 


di<] 


Apply  the  Extended  Hamilton's  Principle:  f(5L  +  5  W)dt  +  BCs  =>  PD£  equations  of  motion 


shown  in  figure  2.1,  and  the  notation  for  the  motion  variables  are  given  in  figure  2.2.  We  seek  to 
cany  out  large  angle  rotational  maneuvers  and  vibration  arrest  with  a  single  control  torque  u(t ) 
input  which  acts  on  body  0  about  the  vertical  axis.  The  resulting  hybrid  system  of  ordinary  and 
partial  differential  equations  of  motion  are 

I  hub  44  =  w  +  E  H,  ,  H,  ={M0  ~S0 10  )i 

at  i=i 


Hi - J  p‘*i  +  X‘  +  m‘  (  A  +  ~^T  \l,  )  +HOT 


d7&  ,  ¥y,i 


P ^ ~d7^ +  ^~dx?  “  0+HOT  •  ‘  =  i.2, 3. ■» 


(2.2.11) 


HOT  indicates  other  known  linear  &  nonlinear  effects  (such  as  rotational  inertia  effects,  rota¬ 
tional  stiffening,  foreshortening  effects,  shear  deformation,  etc.).  The  boundary  conditions  are 


at  x,  =/„:  y,(t,  lj  =  0, 

r)2 

at  x i  -  l. :  h  =  0  +  hot, 

ox, 


(2.2.12) 


The  total  energy  of  the  system  (constant  in  the  absence  of  control  or  disturbances)  is: 

2E  3S  +  &I,)1]  (2.2.13) 


In  view  of  the  above  energy  integral,  we  investigate  the  Liapunov  function 

2U  -  «./»<§)*  *  ♦*$)*«!*  *  /siitfiS)’**  +  ♦£k)*] 

4  *«5(e-0/V 


(2.2.14) 


The  positive  weighting  coefficients  a,  >  0  allow  relative  emphasis  upon  five  substructures’ 
contributions  to  the  total  error  energy  of  the  system.  Note  that  the  open  loop  system  energy 
integral  of  Eq.  (2.2.13)  does  not  depend  upon  the  rigid  body  displacement,  the  final  term  is 
introduced  so  that  Eq.  (22.14)  has  its  global  minimum  at  the  targetfinalstate: 

{ 0.  G}*,.,*,  =  (0/.  0) .  {  yi (*. /),  }da  w  =  ( 0, 0),  f=  1. 2. 3. 4 

More  generally,  error  energy  can  be  measured  from  a  time  varying  target  trajectory  [jl]. 

Differentiation  of  Eq.  (2.2.4),  substitution  of  the  equations  of  motion  (Eqs.  (2.2.1),  (2.2.2)),  and 
substantia]  calculus  leads  to 
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jj  j  4 

V  =  -^-=0  [aoM  +85(0-0/)  +Z(ai-a0)|Xi]  (2.2.15) 

Alternatively,  and  much  more  efficiently,  we  could  have  written  Eq.  (2.2.15)  down  immediately 
using  free  body  diagrams  and  Eq.  (  2.2.8). 


Since  we  require  that  V  <  0,  we  set  the  [  ]  term  to  -  6  and  this  leads  to  U  =  -  &2  and  the 

control  law: 

u  -  — —  [  a5(0-0/)  +a60+  S(at  -ao)|X,  j  (2.2.16) 

3o  i=l 

which,  again,  we  could  have  written  directly  from  Eq.  (2.2.10).  Thus  we  see  that  the  following 
linear,  spatially  discrete  output  feedback  law  satisfies  the  sufficient  condition  ( 0  <  0)  to  globally 
stabilize  this  distributed  parameter  system: 


M  =  -[s,(Mf)  +g2Q+  Sg.-lAi] .  gffi  !>0,  0,  giSh±±>>-  x  (2-2.17) 

2 

The  pervasive  dissipation  condition  that  0  =  -2tf,  6  is  strictly  negative,  for  asymptotic  stability, 


is  satisfied  only  if  the  system  is  fully  controllable.  In  the  linear  case,  we  find  that  the  anti¬ 
symmetric  in  opposition  modes,  (for  a  perfectly  symmetric  structure,  4  identical  appendages) 
have  zero  hub  motion  &  are  uncontrollable  by  a  hub  torque  actuator  —  however,  these  modes  are 
also  theoretically  (assuming  perfectly  identical  appendages  and  clamped  boundary  conditions) 
un-disturbed  for  rest-to-rest  maneuvers  using  a  hub  actuator  [26].  It  is  of  significance  that  we 
have  proven  that  the  same  results  [Eqs.  (2.2.15)  -  (2.2.17)]  are  obtained  when  we  generalize  the 
physical  modeling  assumptions  above  to  include  any/all  of  the  following  effects:  (1)  shear 
deformation  and  rotary  inertia,  (2)  rotational  stiffening  and  foreshortening  effects,  (3)  any/all 
positive  semi-definite  functionals  modeling  the  mechanical  potential  energy  storage  associated 
with  beam  deformation. 


The  invariance  of  the  form  for  the  stabilizing  control  law  and  stable  gain  region,  with  respect  to 
the  most  common  variations  in  modeling  assumptions,  represents  an  important  generalization  of 
the  well  known  robustness  obtained  using  only  the  positive  local  velocity  feedback  term  in  Eq. 
(2.2.17).  The  physical  source  of  this  robustness  is  the  truth  that  the  energy  rate  is  always  given 
by  Eq.  (2.2.6),  irregardless  of  whether  or  not  we  have  correctly  modeled  the  actual  system’s 
physics.  Therefore  we  have  essentially  restricted  the  discussion  to  control  laws  which  cause  the 
error  energy  to  decrease.  The  set  of  system  physical  models  which  will  be  stabilized  at  the  target 
state  are  such  that  the  physically  correct  energy  substructure  functionals  £(,  when  substituted  into 
Eq.  (2.2.1),  yield  a  positive  definite  functional  with  its  global  minimum  at  the  target  state.  Thus 
by  rigorously  stabilizing  a  large  family  of  system  models  by  the  same  control  law,  we  are  then 
left  to  hope  that  the  set  of  stabilized  models  contains  our  actual  system’s  physics.  It  is  fortuitous 
that  almost  all  beam  constituitive  modeling  assumptions  indeed  result  in  a  positive  beam  poten¬ 
tial  energy  functional  which  vanishes  for  zero  deformation,  and  obviously,  zero  deformation  is 
the  most  common  target  state!  Of  course  the  accuracy  of  the  predicted  closed  loop  response 
(while  stability  may  be  assured)  is  indeed  a  function  of  the  degree  to  which  the  system  is  accu¬ 
rately  modeled,  and  clearly  the  optimization  of  the  control  gains  is  model-dependent. 
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For  the  special  case  for  which  the  four  appendages  and  tip  masses  are  assumed  identical,  we  have 
studied  the  above  control  law  and  the  generalizations  thereof  analytically,  numerically,  and  have 
conducted  successful  laboratory  experiments,  [see  reference  26  and  Attachment  3] 


2.2.3  Status  and  Outlook 

As  is  evident  from  the  above  discussion  and  Attachment  1,  we  have  made  significant  progress  on 
several  fronts:  (1)  development  of  a  novel  approach  for  designing  stable  control  control  laws  for 
nonlinear  distributed  parameter  systems,  (2)  analytical  and  numerical  studies  of  the  validity  of 
the  approach  for  special  cases,  and  (3)  laboratory  experimental  validation  of  the  approach  for  a 
simple  multi-body  maneuver  configuration.  The  results  to  date  have  been  very  encouraging.  We 
have  also  addressed  extensions  of  the  basic  methodology  to  permit  tradeoffs  between  competing 
measures  of  optimality  such  as  minimum  maneuver  time  versus  minimum  vibration  measures 
(see  attachment  3).  During  the  next  year,  we  expect  to  devote  most  of  our  effort  to  three  theoreti¬ 
cal  issues:  (i)  Addressing  the  issues  raised  by  uncontrollable  and/or  poorly  controllable  dynami¬ 
cal  sub-spaces  of  motion  (e.  g.,  for  a  linear  structure,  uncontrollable  natural  vibration  modes)  for 
nonlinear  systems  as  a  function  of  actuator  configurations,  (ii)  Extending/modifying  the  power 
principle  to  consider  quasi  coordinate  descriptions  of  the  system  kinematics,  (iii)  Extending/ 
modifying  the  power  principle  to  address  non-natural  systems  for  which  the  kinetic  energy  is  not 
a  symmetric  quadratic  form  in  the  generalized  velocities,  and  the  potential  energy  is  not  a  posi¬ 
tive  definite  function  of  the  generalized  coordinates.  All  three  of  these  theoretical  issues  will  be 
studied  in  the  light  of  simple  example  structures,  motivated  by  potential  practical  applications, 
with  appropriate  characteristics  to  illuminate  the  salient  features  of  the  developments. 
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2.3  Symbolic  Derivation  of  Nonlinear  Feedback  Control  Laws 

2.3.1  Basic  Ideas 


In  a  somewhat  unorthodox,  but  we  trust,  effective  format,  we  motivate  the  symbolic  nonlinear 
control  design  approach  by  first  presenting  several  example  nonlinear  control  problems  then  their 
solutions  which  we  have  recently  been  able  to  carry  to  completion.  We  first  state  the  three 
problems,  then  outline  the  main  features  of  their  solutions.  The  first  two  problems  are  overly 
"academic"  examples  selected  because  they  offer  a  transparent  way  to  introduce  these  ideas.  The 
third  problem  has  obvious  practical  significance,  this  example  illustrates  the  use  of  this  approach 
to  design  an  optimal  feedback  control  law  for  the  nonlinear  spacecraft  attitude  maneuver 
problem.  These  ideas  grew  from  our  historical  research  documented  in  references  [29,  30]. 

Three  Problem  Statements 


To  illustrate  the  concepts,  we  introduce  the  differential  equations  and  performance  indices  for 
three  nonlinear  control  problems.  The  first  two  problems  are  provided  for  a  simple  introduction, 
but  they  have  a  similar  structure  to  the  third  problem  (large  angle,  nonlinear  spacecraft  attitude 
maneuvers),  and  the  same  methodology  readily  solves  all  three  problems. 


Problem  I 

Find  an  optimal  feedback  control  law  to  compute  u(x)  for  the  system  described  by 
x  =  -x  +  ax2  +  u  ,  a  >  0 
which  minimizes  the  performance  measure 

/  =  i  \  (x2  +  u2)  dt 
z  o 


(2.2.18) 


(2.2.19) 


Problem  II 

Find  an  optimal  feedback  control  law  to  compute  Uj(xjf  x2),  u2(x},  x2)  for  the  system 
Xi  =  -  X]  +  X]X2  +x%  +  Ui 

2  (2.2.20) 

X2  =  -  x2  +  X]X2  +  X]  +  u2 


1  rf  ^ 

which  minimizes  the  performance  measure  J  =  ^  J  X  ( +  u2)  dt  (2.2.21) 

Z  0  i- 1 

Note  that  Eqs.  (2.2.20)  can  be  written  in  an  alternate  matrix  format  as 

Jc=AiXi  +A2X2  +Bu  ,  (2.2.22) 

with  Xrl  =  xr=[xl  x2],x[=[x2i  X]x2 
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As  we  show  below,  the  above  notation  generalizes  to  encompass  high  dimensional  nonlinear 
systems.  The  linear  appearance  of  the  above  notation  should  not  obscure  the  fact  that  the  system 
of  Eq.  (2.2.22)  is  nonlinear,  note  that  the  vector  x2  contains  all  of  the  quadratic  nonlinear  terms. 


Problem  III 

Find  an  optimal  feedback  law  for  the  control  torques  u{qv  q2,  q3,  co;,  to2,  0)3),  for  i  =  1 ,2,3  to 
control  the  six-state,  three-input  system  described  by  [  30]  the  following  system  of  equations 

<7i  =2  [(l+<7?)G>i  +(<7i  <72 -<73)0)2  +  (<7i  <73 +<72  )&)3  ] 

<72  -  5  [(<71  <72  +<?3  )C0i  +  ( 1  +<72  )o)2  +  (<72  <73  ~<7i  )c»3  ] 

<73  =  ^3~<?2)0)l  +  (^2<73+<7l  )t»2  +(1 +<73)0)3] 

2  (2.2.23) 

0)1  =(V1)(°2(1>3  +(i)«l 

6)2  =(V)(°3a)'  +(^)M2 

(i>3  =(L^)0)I0)2  +(^)«3 


which  minimizes  the  performance  index  /  {.X1  Qx  ~^~UTRu)dt  (2.2.24) 

0 

0,  and  /?  are  symmetric  positive  definite  weight  matrices;  the  state  and  control  vectors  are 
X  =  [ <71  <72  <?3  0)1  0)2  (03],«r  =  [Ml  U2  M3] 

The  q’ s  are  the  three  Rodriguez  attitude  coordinates  [28]  and  the  co’s  are  the  three  orthogonal 
components  of  the  angular  velocity  vector  along  principal  axes.  The  u’s  are  the  control  torques 
about  the  principal  axes.  The  /’ s  are  the  principal  inertias.  These  equations  govern  the  general, 
large  angle,  nonlinear  attitude  maneuvers  of  a  rigid  spacecraft.  An  advantage  of  *he  Rodriguez 
parameters,  compared  to  any  choice  of  three  Euler  angles,  is  that  no  transcendental  functions 
appear.  The  above  equations  are  exact,  the  degree  of  polynomial  nonlinearity  is  three.  Note  that 
the  target  orientation  of  the  body  (qt  =  0)  has  been  selected  as  the  inertial  frame. 

Equations  (6)  can  be  viewed  as  a  special  case  of  the  most  general  nonlinear  differential  equation 

J t  =  ZAiXi  +  Bu  (2.2.25) 

1=1 

where  xi  are  column  vectors  containing  the  6  linear  state  variables  (x  =  x,),  the  21  distinct 
quadratic  combinations  of  the  state  variables  (x2),  and  the  56  cubic  combinations  of  the  state 
variables  (x3),  for  above  case  (as  in  all  others  studied  to  date!),  the  three  Ai  matrices  are  very 
sparse,  as  given  in  detail  below. 
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A 

*1*2 
*1*3 
JCi  JC4 
*1*5 
*1*6 

A 


r*r 

*2 


*2*3 

*2*4 

*2*5 


I  2111 

> ,  *2  =  *2*6 


*4 

*5 

.*6> 


*3*4 

*3*5 

*3*6 


*4 


56(1 

*3 


*4*5 

*4*6 


A 


*5*6 

4 


■i 

*1*2 
A  *3 
*?* 
*i*5 
*?*« 
*1*2 
*1*2*3 
*1  *2^ 
*1*2*5 
*1*2*6 
*1*5 
*1*3*4 
*1*3*5 
*1*3* 
*1*4 
*1*4*5 
*1*4*6 
*1*5 
*1*5*6 
*1*6 
4 

4* 

*2*4 

*2*5 

*2*6 

*2*3 

*2*3*4 

*2*3*5’ 

*2*3*6 

*2*2 

*2*4*5 

*2*4*6 

*2*2 

*2*5*6 

2 


*5*1 

4*6 

*»4 

*l*«*6 

•>*5 

*1*5*6 

*j4 

4 

*5*1 

4*6 

*•*5 

*1*5*6 

*»4 

4 

4*6 

*s4 

4 


The  non-zero  elements  of  the  A,  are: 
Aj(l,4)=A,  (24)  =  A,  (3,6)  =  | 


6(21  6k2I  6x21  1 

A2(l,ll)  =  A2(2.13)  =  A2(34)  =  2 

6i2I  6x21  6x21  | 

A  2  (3,14)  =  A  2  (2,6)  =  A  2  (3,9)  =  — ^ 

A  2  (4,20)  =  (n^1) 

A2(5,18)  =  (^) 

A2(6,17)=(^) 


A3(l,4)  =  A3  (1,10)  =  A  3  (1,15)  = 


fa56  4b56  6r56 

A  3  (2,9)  =  A  3  (2,25)  =  A  3  (2,30)  = 


6x56  6c  56  6i56 

A  3  (3,13)  =  A  3  (3,29)  =  A  3  (3,37)  = 
The  B  matrix  is: 


1 

T 

1 

5 

1 

T 


B  = 


0 


0 


0  0 

h  0 
0  i 
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(2.2.26) 


Prior  to  presenting  the  solutions  for  the  feedback  controls  for  the  above  three  problems,  we  state 
that  we  have  developed  a  general  approach  which  applies  the  optimal  control  necessary  condi¬ 
tions  (Pontryagin’s  Principle  and  Pontryagin’s  necessary  conditions  [28]),  to  minimize  an  index 
of  the  form  of  Eq.  (2.2.24))  and  generates  symbolically  the  necessary  conditions,  and  solves  these 
for  the  particular  set  of  symbolic  differential  and  algebraic  equations  governing  the  optimal 
nonlinear  feedback  control  gains  for  any  particular  member  of  the  class  of  nonlinear  dynamical 
systems  described  by  differential  equations  of  the  form 

M 

£ =ZAiXi+Bu  (2.2.27) 

i=i 

where  x  is  an  nxl  state  vector,  u  is  an  mxl  control  vector,  M  is  the  highest  degree  of  polynomial 
nonlinearity,  xt  is  a  column  vector  containing  all  distinct  polynomial  combinations  (of  degree  i) 
of  the  elements  of  the  state  vector  x.  It  is  obvious  that  Eqs.  (2..2.27)  and  the  generalized  problem 
statement  embraces  a  very  large  class  of  systems,  including  all  three  of  the  above-stated  prob¬ 
lems  as  special  cases,  and  therefore  a  large  family  of  nonlinear  mechanics  problems  arising  in 
structural  dynamics  and  control.  The  nonlinear  feedback  controls  u  are  found  as  direct  solutions 
of  the  general  necessary  conditions  for  optimal  control  of  the  general  family  of  systems  described 
by  Eqs.  (2.2.27).  Our  approach  leads  directly  to  symbolic  equations  satisfied  by  the  optimal 
feedback  gain  matrices  Gt  in  the  polynomial  expansion 

N  N 

U  =  "LG iXi ,  where  G,  =  -  R'1  BT Ki,  X  =  (2.2.28) 

1=1  i=l 

The  degree  (AO  of  the  nonlinear  feedback  control  expansion  is  not  restricted  to  be  equal  to  the 
degree  (M)  of  the  nonlinearity  in  the  original  differential  equations  (it  is  fortuitous,  as  is  con¬ 
firmed  in  the  examples  below,  that  N  typically  required  for  practical  convergence  is  in  fact 
usually  of  low  degree,  but  no  universal  conclusion  can  be  made).  Note  that  X  is  an  nxl  vector  of 
Lagrange  multipliers  which  arise  in  the  optimal  control  necessary  conditions  [30]. 

It  is  significant  to  note  that  we  can  find,  via  symbol  manipulation,  a  set  of  sequentially  solvable, 
general  algebraic  and/or  differential  equations  satisfied  by  the  gains  K{  where  the  matrices  Ai?  B, 

Q,  R  ,  or  subsets  thereof,  can  appear  as  algebraic  parameters  (i.  e.,  we  do  not  first  have  to  first 
specify  numerical  values  for  the  system  parameters,  or  even  numerical  elements  for  the  weight 
matrices);  we  can  leave  any  subset  (or  all)  of  the  system  parameters  as  symbols  and  develop  a 
sequence  of  equations  which  are  specialized  for  the  size  and  sparsity  patterns  of  the  matrices  of 
the  particular  system  of  interest.  These  equations  can  then  be  solved  numerically  for  the  gains 
corresponding  to  the  particular  applications  of  interest.  Another  way  of  viewing  our  result  is  to 
observe  that  it  has  been  known  for  several  decades  that  the  optimal  linear  control  gains,  for  a 
linear  system,  are  generated  through  solution  of  a  matrix  Riccati  equation  which  depends  ex¬ 
plicitly  upon  matrices  A,  B,  Q,  R  ;  we  have  conceived  of  a  systematic  way  to  develop  the 
analogous,  explicit,  sequentially  solvable  equations  governing  the  higher  order  nonlinear  feed¬ 
back  gains,  for  the  class  of  systems  described  by  Eq.  (12.2.27). 

It  is  apparent  that  a  general  realization  of  this  approach  would  be  very  attractive,  not  only  be¬ 
cause  it  makes  determination  of  the  equations  governing  the  nonlinear  gains  relatively  routine, 
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but  also  because  changes  in  subsets  of  the  system  parameters  or  weights  in  the  performance 
measure  can  be  quickly  accounted  for  by  simply  changing  the  desired  parameters  in  the  (one¬ 
time-derived)  gain  equations  and  then  simply  re-solving  these  equations  for  the  new  gains.  This 
approach  would  be  quite  impossible,  for  nonlinear  systems  of  even  moderate  dimensionality, 
without  use  of  modem  computer  algebraic  manipulation.  While  a  given  dynamical  system  might 
submit  after  man-months  of  algebra  and  calculus,  it  would  only  rarely  be  implemented  success¬ 
fully  due  to  the  associated  large  investment  of  human  effort,  elapsed  time,  and  the  intimidating 
problems  raised  by  debugging  the  results  and  integrating  this  process  into  an  invariably  iterative 
controller  design  cycle.  It  appears  evident  that  a  successful  implementation  of  our  approach  will 
make  possible  routine  application  of  perturbation  methods  to  derive  nonlinear  control  laws,  at 
least  for  large  families  of  problems  (not  merely  in  principle ,  but  in  fact\). 

Our  first  implementation  (  using  macsyma  )  of  these  ideas  can  be  applied  "fairly  routinely" 
without  requiring  intolerable  storage  or  execution  time,  if  the  product  of  MN  is  less  than  about 
20.  Due  to  the  several  curses  of  dimensionality,  we  do  not  feel  that  this  approach  will  prove 
practical  in  the  near  term  if  MN  is  greater  than  about  40,  but  this  is  still  encompasses  a  very  large 
class  of  problems.  In  the  applications  to  date,  we  have  found  the  matrices  involved  and  the 
resulting  control  gains  are  rather  sparse,  it  is  possible  that  the  curse  of  dimensionality  can  be 
substantially  reduced  by  introducing  (as  yet  undeveloped)  methods  to  anticipate  and  take  advan¬ 
tage  of  the  sparse  structure  of  the  particular  problem.  Immediate  extensions  to  higher  dimensions 
can  be  undertaken,  but  we  feel  that  effects  and  problems  associated  with  dimensionality  increases 
should  be  studied  carefully  in  the  context  of  a  systematic,  escalator-styled  research  effort,  taking 
time  to  consider  specific  applications.  This  will  permit  the  evolving  formulations  and  computer 
implementations  to  benefit  fully  from  the  insights  which  stem  from  using  the  methods  on  prob¬ 
lems  small  enough  to  make  the  salient  features  transparent.  Since  this  approach  represents  a  new 
controller  design  methodology,  the  basic  research  will  no  doubt  benefit  from  the  analytical  and 
artistic  insights  gained  from  several  case  study  applications.  Convergence  proofs  are  not  avail¬ 
able  for  arbitrary  systems  belonging  to  the  general  class  described  by  Eqs.  (2.2.27  -  2.2.28),  but 
convergence  will  be  studied  on  a  case-by-case  basis,  and  we  will  seek  to  establish  insights  on 
how  to  approach  resolving  the  convergence  issue  for  high-dimensioned,  high-order  expansions. 

In  order  to  make  the  essential  ideas  easy  to  understand  and  to  simultaneously  display  some  of  the 
progress  we  have  made  to  date,  we  present  the  solution  details  for  the  first  problem  as  we  origi¬ 
nally  developed  them  by  hand ,  these  same  results  have  been  confirmed  and  extended  using  the 
symbolic  manipulator.  For  the  second  and  third  problems,  we  provide  only  solution  outlines  and 
some  of  the  ’end  products’  of  this  development,  due  to  space  limitations.  All  three  of  the  above 
stated  problems  have  been  solved  by  hand  (to  low  order)  to  verify  the  correctness  of  the  symbol 
manipulation  implementation,  we  have  also  successfully  completed  comparisons  of  the  results 
with  Carrington’s  dissertation  [29].  We  now  discuss  the  solutions  of  the  three  problems. 
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2.2  Solution  of  Problem  I 

We  seek  an  optimal  feedback  control  law  to  compute  u(x)  for 

x  =  -x  +  C XX1  +u,  a  >0  (2.2.29) 


which  minimizes  the  performance  measure  J  ~  \ 

The  Hamiltonian  for  this  system  is 

H=  ^  (x2  +  u2 )  +  X(-  x  +  a  x2  +  u) 


J  (jc2  +  u2)  dt 
o 


(2.2.30) 

(2.2.31) 


The  Pontryagin  necessary  conditions  are: 


rdH 


=  0 


u  =  -  X 


< 


X  =  - 


dH 


~f  A,  —  2 qlx\  >(2.2.32) 


_  ZH 
x  -  IK 


=  -x  -  X  +  a  x2 


J 


We  seek  a  feedback  form  of  the  control  law,  specifically,  we  seek  the  polynomial  gains  Ki  in 
N 

-u=X='LKlx‘  (2.2.33) 

i  =  1 


Substituting  Eq.  (2.2.33)  and  its  time  derivative  into  Eqs.  (2.2.32),  we  find  the  homogeneous 
condition 

[£,  -  2KX  -  k\  +  \]x  +  [k2  -3(1  +Ki  )K2  +30A:,  ]x*  + 

[*s  -4(l+tf,)*3  -+-4ocATj  -2k\]>?  +  ...  +  [kN  ~(N+l){l+kx)kN  -FN(aJ<x , ....  kN-x  )]x"  =0 


Since  the  above  equation  must  hold  at  every  point  in  the  state  space  (i.  e. ,  for  all  x),  we  conclude 
that  all  [  ]’ed  coefficients  must  vanish  independently,  this  provides  the  following  equations: 


ki  -  2k x  -  k]  +  1  =  0 

kx 

k2  -3(1+*1  )k2  =  -  3<xATi 

=> 

k2 

*3  -4(l+fl(i  )k2  =  -  AaK\  +  2 k\ 

k2 

(2.2.34) 

kN  -(N+lKl+AT,  )kN  =  Fn( a,  k},k2, ... X n-i ) 

=> 

ks 

Notice  the  following  structure  and  properties  of  Eqs.  (2.2.34)  satisfied  by  the  optimal  gains: 
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®  These  equations  can  be  solved  sequentially  for  the  A'-,  all  but  the  first  are  linear  equations. 

®  The  equations  can  solved  to  arbitrary  order,  since  we  have  developed  explicit  recursions 

for  Fn. 

®  The  first  equation  for  K,,  is  a  scalar  Riccati  equation  (no  surprise  here!). 

®  If  we  impose  a  =  0  and  the  boundary  conditions  Rfo)  =  0  [consistent  with  X.(tf)  =0  as  a 

transversality  condition  for  free],  we  see  that  all  nonlinear  gains  vanish  identically  and 
therefore  the  above  is  indeed  a  direct  generalization  of  the  classical  linear  regulator  with  a 
quadratic  performance  index,  the  optimal  control  of  Eq.  (2.2.33)  reduces  to  u  =  -  x. 

®  If  tf  — >  ,  analogous  to  the  classical  steady  state  regulator,  we  can  show  that  all  Ki  ap¬ 

proach  constants  and  Eqs.  (2.2.34)  therefore  reduce  to  a  sequence  of  algebraic  equations 
for  the  constant  gains  (we  can  show  that  the  positive  real  root  is  the  proper  selection  for 
Kj),  and  the  solution  for  the  higher  gains  involves  only  simple  algebraic  operations.  For 
this  case,  with  a  =  1,  the  first  few  gains  are  numerically: 

K,  =  0.41321,  K2  =  0.29289,  ^  =  0.17678,  AT,  =  0.088388,  =  0.033145. 

Shown  below  in  Figure  2.3  is  the  performance  index  versus  the  order  N  of  the  feedback  control 
and  the  trajectories  of  jc(t)  and  u(t)  for  typical  initial  conditions  [jc(0)=1.3].  As  is  qualitatively 
evident,  the  nonlinear  terms  are  constructive  and  convergence  is  rapid.  We  have  shown  that  the 
nonlinear  controls  for  N>3  are  globally  stable,  whereas  large  but  finite  domains  of  stability  are 
associated  with  the  linear  and  quadratic  feedback  control  laws.  This  solution  has  been  verified 
by  the  symbolic  manipulator  code.  As  we  indicate  below,  the  structure  of  this  simplest  scalar 
example  fully  generalizes,  so  that  the  above  observations  apply  in  a  much  more  general  context. 
However,  the  stability  and  convergence  issues  are  problem-dependent,  as  should  be  expected. 
We  have  established  a  means  to  generate  the  matrix  equivalent  of  the  sequence  of  scalar  equa¬ 
tions  ( Eqs.  (2.2.34  )  and  their  solutions  for  the  optimal  nonlinear  feedback  gains 

Solution  of  Problem  II 

We  seek  an  optimal  feedback  control  law  to  compute  u/Xj,  x2),  u2(xv  x2)  for  the  system: 

X  i  =  -X]  +X/X2  +X%  +  U] 

(2.2.35) 

X-2  =  -  X2  +  XjX 2  +xj  +  U2 
which  minimizes 


J  =  j  +«!)^  =  2 )(xrQx  +  utR u)dt 

0  *“1  n 


(2.2.36) 


Eqs.  (2.2.35)  can  be  written  as  Jc  =  A\X\  +  A2X2  +  Bu  , 


Q=R  -1 


(2.2.37) 
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where  x[=xT=[*i  x2],x\  =  *1*2  xf],ttT=[«,  «2],Ai  =  [q  °] ,  ^2  =  [g  j  i],fi  =  [o  1]' 

The  Hamiltonian  function  is  H  =  ^  (x7^ JC  +  uT u)  +  XT  (A\Xj  +A2X2  +  Bu)  (2.2.38) 

dH 


X  =  - 


The  Pontryagin  conditions  are  [5]  J 


Hu 

dH 

Hx 

dH 


=  0 


u  =  -BTX 


\  (2.2.39) 


x  —  —  A 1  Jiff  +  A2x2  —  B  B  X 

—[S]  ■[:::] 


J 


We  seek  to  determine  the  feedback  gains  Kt  in  the  control  law  expansion 

u=-BT\=!LGixi,  where  X  =  ZKiXi  ,  G,  =  -  BT  Kt  _ 

1=1  J=1  (Z.Z.4U) 

Substitution  of  the  expansion  of  Eqs.  (23)  into  the  necessary  conditions  of  Eq.  (24)  and  collecting 
terms  leads  directly  to  a  system  of  two  homogeneous  conditions  of  the  following  form 


FuiKxyx,  +  Fa  (Kl]x2+Fni(K\,K2'rf  +  Fni  (AT, ,  K2  ]x2  jc2  +  Fm(Ki,  K2)£  +  ...  =  0,/  =  l,2 

Requiring  these  conditions  to  hold  at  every  point  in  the  state  space,  we  can  set  the  coefficients  of 
all  powers  and  products  of  the  elements  of  x  to  zero.  Upon  carrying  through  the  algebra,  we  find 
the  linear  terms  yield  four  algebraic  equations  Fu  (K\ )  =  F\2(K\ )  =F2i  )  =  F22(Ki )  =  0  , 
which  are  precisely  the  four  elements  of  the  Riccati  equation  [29,  30] 

K\A\  +  A\Kx  - K\BR-' Bt K\  +Q  =0  (2.2.41) 

The  Riccati  equation  can  be  solved  for  the  symmetric  linear  gain  matrix  Kx.  Setting  the  six  quad¬ 
ratic  term’s  coefficients  to  zero  yields  FUl (Kx ,  K2)  =  Fx2t  (K\  ,K2)=  F22i (ATj ,  K2 )  =  0  ,  i  =  1,  2; 
we  find  that  these  six  algebraic  equations  are  linear  in  the  six  distinct  elements  of  K2 ,  and  can  be 
brought  to  the  form  of  the  linear  system 

[IjtiKx))  vec{^2>  =  R2 (K\ )  (2.2.42) 

where  is  a  6x6  matrix  whose  elements  are  functions  of  K\,  R2 (ATj )  is  a  6x1  vector 

whose  elements  depend  upon  K\,  and  vec{K2 }  is  a  6x1  vector  whose  elements  are  the  six  distinct 
elements  of  K2.  Obviously,  if  [L2  (K\ )  ]  is  of  full  rank,  Eq.  (2.2.42)  can  be  inverted  for  vec{/f2 }. 
To  conserve  space,  we  do  not  write  out  the  algebraic  equations  for  the  elements  of  these 
matrices.  We  find  that  Eq.  (2.2.42)  generalizes;  the  higher  order  gains  are  determined  by  a 
sequence  of  linear  equations  of  the  form 
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IU(KX .  K2, ....  Kk-X)\  vec{Ari}  =  R*(ATi ,  K2, Kk.x ),  k  =  2,3,...  N 


(2.2.43) 


For  the  case  of  particular  Ait  B,  of  Eqs.  (2.2.20),  and  Q=  R  =  /,  we  have  carried  through  the 
above  developments  and  find  the  following  numerical  values  for  the  first  five  control  gains  (G,= 
-R%): 

q  __r0.41421  0  1  „  _  _r  0  0.39052  0.195261  _  __ TO. 12459  0.27022  0.22222  0.09007 T 

1  L  0  0.41421  J’  2  L0.19526  0.39052  0  J  ’  °3  [.0.09007  0.22222  0.27022  0.12459  J 


g  =_r°05125  017517  026213  017475  oo4379i  c  =  _r°°1656  008537  016605  016225  008302  °oi7o7i 

UA  [.0.04379  0.17475  0.26213  0.17517  0.05125J’  US  L°  01707  0.08302  0.16225  0.16605  0.08537  0.01656 J 

Since  this  system  is  highly  nonlinear,  simply  ignoring  the  nonlinear  terms  and  deriving  an 
approximate  linear  optimal  control  law  will  be  valid  only  near  the  origin.  In  Figure  2.4  we  show 
graphs  of  the  stable  region  (shaded)  vs  N  and  typical  trajectories  of  the  state  and  control  variables 
for  typical  initial  conditions,  for  controllers  based  upon  linear  (N=  1)  through  quintic  (N=5) 
feedback.  Notice  (as  might  be  anticipated),  the  linear  feedback  law  does  not  globally  stabilize 
this  nonlinear  system,  but  it  does  near  the  origin  of  the  state  space.  Including  the  quadratic  terms 
modifies  the  stable  region,  increasing  the  stable  domain  area  in  the  positive  (xj ,  x2 )  quadrant, 
but  decreasing  it  elsewhere.  Including  the  third  degree  terms  results  in  global  stability. 

There  is  evidence  that  outside  the  shaded  region  of  Figure  2.4b,  all  even  degree  feedback  control¬ 
lers  are  unstable,  whereas  inside  this  region  all  controllers  with  2<N<5  are  stable  and  converge 
rapidly  to  the  optimal  control.  Using  Lyapunov  methods,  we  determined  an  N=2  globally 
stabilizing,  sub-optimal  quadratic  feedback  law  for  this  system,  but  this  does  not  detract  from  the 
significance  of  the  above  results,  since  they  generalize  fully  to  the  case  of  higher  order  systems 
with  a  small  number  of  controllers,  in  which  case  no  general  method  exists  to  construct  a  control 
law  for  which  global  stability  is  guaranteed.  Indeed,  we  believe  the  ideas  we  are  discussing  can 
be  developed  into  a  widely  applicable  method  for  enlarging  the  stable,  near-optimally  controlled 
region  for  a  large  class  of  nonlinear  systems.  It  is  evident  that  including  the  nonlinear  terms  in 
the  truncated  feedback  control  law  is  constructive  and  convergence  to  near  optimal  solutions  can 
be  achieved  over  greatly  enlarged  regions  of  the  state  space.  However,  careful  evaluation  of  the 
solution  behavior  is  still  required,  there  are  no  guarantees  of  monotonic  convergence. 
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Figure  2.4  Stability  Maps  and  Typical  T rajectories  for  Problem  II: 

Effects  of  Increasing  the  Degree  of  the  Nonlinear  Feedback 


Solution  of  Problem  III  (Optimal  Nonlinear  Feedback  Control  of  Spacecraft  Attitude) 

We  seek  an  optimal  feedback  control  law  to  compute  the  feedback  torque  u(x)  to  maneuver  the 
system  modeled  by  Eqs.  (2.2.23)  or  more  generally  Eqs.  (2.2.25)  to  the  origin  of  the  state  space 
(x  =  0)  in  such  a  fashion  that  minimizes  the  index  of  Eq.  (2.2.24).  The  detailed  discussion  of  this 
example  is  not  given,  but  is  analogous  to  the  solution  of  Problem  II.  The  expansions  were 
carried  out  to  3rd  order  via  the  symbolic  manipulator.  This  is  essentially  the  same  system 
considered  in  [29,30].  Here  we  present  only  a  summary  of  the  numerical  solution  of  the  resulting 
general  symbolic  gain  equations  for  the  particular  system  parameters:  /7  =  1.00,  /2  =  0.83,  I3  = 
0.92,  Q  ,  R  identity  matrices,  and  unit  initial  conditions  on  all  state  variables. 


For  the  First  two  gain  matrices  (G,  =  -  R1BrKi  ),  where  the  AT-  are  from  the  expansion  of  Eq. 
(2.2.28),  we  find  the  following  numerical  values 
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Due  to  space  limitations,  we  show  only  a  few  digits  and  do  not  display  the  third  order  gain 
matrix  which  was  also  computed.  As  is  apparent,  the  particular  system  dynamics,  inertias,  and 
identity  weight  matrices  resulted  in  sparse  gain  matrices,  this  pattern  carries  over  to  the  cubic 
gains.  This  sparsity  pattern  is  unaffected  by  variations  of  the  elements  in  the  diagonal  inertia,  Q, 
and  R  matrices,  although,  obviously,  the  numerical  values  of  the  gains  are  affected.  The  sparsity 
pattern  of  the  optimal  gain  matrices  allows  us  to  write  down  the  form  of  the  feedback  control  law 
explicitly,  as  follows: 

Ui  = 

u2  = 
u3  = 

linear  feedback  terms  quadratic  (note  "gyroscopic"  structure)  feedback  terms  3rd  &  HOT 


to,  (l.i)qi  +  g,  (i,4)b)i :  + 
G,(2^>q2  +  G,  (2,5)02:  + 
G,(3.3)q3  +  G,  (3,6)03;  + 


bVoijqi  q3"  +  Gj(i,i  i) q2  03'  +  'oYa.i-bqYaY  +  cfc(i‘jo)Ct^  ojj 

G,(2,3)q3qi  +  Gj (2,13)q3  CO]  +  G,(2,6)qi03  +  0,(2,18)030] 
G,(3^)q]q2  +  G,(3,5)qi02  +  G,(3,9)q2Oi  +  G,C3,17)0i02 


:+ 
4- . 


While  it  is  easy  to  conjecture  the  uncoupled  structure  of  the  linear  terms  in  the  optimal  control 
(the  linearized  dynamical  equations  are  uncoupled),  and  after-the-fact,  the  "gyroscopic"  structure 
of  the  optimal  quadratic  feedback  is  intuitively  reasonable,  it  is  more  difficult  to  anticipate  the 
structure  of  the  3rd  and  higher  order  terms  of  the  optimal  control  law.  The  third  degree  feedback 
gains  are  similarly  sparse  (only  15  of  the  56  elements  in  each  row  of  G3  are  non-zero),  but  are  not 
displayed,  for  brevity.  The  beauty  of  the  above  method  is  that  not  only  the  structure  of  the 
control  law,  but  also  the  optimal  values  of  the  gains  can  be  routinely  determined.  Of  course, 
fully  populated  weight  and  inertia  matrices  result  in  more  densely  populated  gain  matrices. 
Shown  in  Figure  2.5  are  graphs  of  the  the  state  and  control  variable  trajectories  corresponding  to 
unit  initial  conditions  (note  the  initial  conditions  correspond  to  large  angular  velocities  about 
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Figure  2.5  Perturbation  Feedback  Controlled  Spacecraft  Attitude  Maneuver 
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each  axis  —  thus  the  gyroscopic  nonlinear  effects  are  fairly  pronounced).  As  is  evident  in  Figure 
2.5,  we  again  find  very  constructive  effects,  upon  augmenting  the  classical  linear  control  gains 
with  the  optimal  quadratic  and  cubic  feedback,  it  is  also  obvious  that  rapid  convergence  is  being 
achieved.  Even  though  the  quadratic  and  cubic  control  contributions  are  significant,  we  antici¬ 
pate  the  fourth  and  higher  order  feedback  contributions  would  make  negligible  contributions. 

These  results  and  those  obtained  in  solving  the  first  two  examples  appear  to  be  typical,  the 
absence  of  formal  convergence  proofs  does  not  usually  prevent  us  from  obtaining  practical 
solutions  which  do  in  fact  display  convincing  evidence  of  convergence  with  low  degree  nonlinear 
feedback.  Especially  important  is  the  promise  of  this  approach  to  make  the  determination  and 
revision  of  nonlinear,  near-optimal  controls  relatively  straight  forward,  so  that  convergence  can 
be  studied  and  we  can  efficiently  incorporate  system  and  performance  index  modifications. 

2.3.3  Status  and  Outlook 

As  is  evident  from  the  above  developments,  we  have  made  some  significant  progress  in  the 
development  and  preliminary  evaluation  of  a  new  approach  to  nonlinear  feedback  control.  We 
have  however  encountered  several  dimensionality  -  related  difficulties.  The  number  of  symbolic 
operations  depends  upon  the  second  power  of  the  product  MN  (M  =  highest  degree  of  polyno¬ 
mial  nonlinearity,  N  =  order  of  the  dynamical  system);  we  have  been  successful  in  carrying  to 
completion  the  design  for  several  examples  for  which  MN  <  20.  However,  it  presently  appears 
that  our  present  formulation  and  implementation  will  not  be  practical  if  MN  exceeds  about  40 
due  to  excessive  computational  demands.  However,  we  are  optimistic  that  new  developments 
which  adaptively  exploit  sparsity  structure  associated  with  each  system  (instead  of  the  present 
most  general  approach  which  initially  allows  for  all  nonlinear  terms  of  degree  M  and  lower)  will 
be  developed  which  will  permit  higher  dimensioned  applications. 

A  more  serious  difficulty,  associated  with  the  lack  of  stability  guarantees,  has  been  encountered. 
We  have  found  several  examples  in  which  nonlinear  controllers  of  certain  degrees  were  unstable; 
this  was  found  after  the  fact  (after  the  optimality  conditions  were  applied  and  the  gains  com¬ 
puted).  Unlike  linear  systems  (for  a  linear  controllable  system,  it  can  be  proven  that  minimizing 
a  quadratic  index  always  leads  to  a  stable  controller),  there  is  no  proof  that  polynomial  (truncated 
at  some  degree)  feedback,  determined  to  minimize  some  positive  performance  index  will  stabi¬ 
lize  a  given  nonlinear  system.  Indeed  we  have  encountered  several  counter-examples.  These 
issues  will  be  studied  further  in  the  next  year. 
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3.0  Concluding  Remarks 

Section  2  summarizes  progress  we  have  made  on  three  sets  of  research  problems: 

Penalty  methods  for  simulation  of  flexible  multibody  dynamics. 

A  power  method  for  design  of  output  feedback  controllers  for  nonlinear  distributed 
parameter  systems 

A  symbolic  approach  for  designing  full  state  feedback  control  laws  for  dynamical  systems 
with  polynomial  nonlinearities,  polynomial 

On  all  three  sets  of  problems,  we  have  obtained  some  fundamental  analytical  results  and  have 
studied  prototype  applications  to  evaluate  salient  features  and  the  practical  potential  of  the 
methodology.  In  all  three  areas  substantial  progress  has  been  made  as  reported  above,  and  as 
discussed  in  the  sub-sections  2.1.4,  2.2.3,  and  2.3.3,  we  are  aggressively  pursuing  extensions  to 
the  results  in  all  three  areas  of  investigation. 
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eneraiization  ot  the  Methodology 

A  Power  Method  for  designing  stable  feedback  controllers 
of  nonlinear,  multi-body  systems  modeled  by  ODE/PDEs. 


Fondest  Hopes  and  Wildest  Dreams:  A  Wish  List 
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The  method  can  accommodate  tradeoffs  between  competing  measures  of 
optimality ,  for  example: 

minimum  time  vs  minimum  control  energy  vs  small  performance  errors 


E  BODY  CONFIGURATION 


p 


i 

i 


j 


z 

o 

H 

o 

§ 

o 

c» 

z 

o 

H 

<C 

D 

O' 

Ed 


P 

3 

.2 

>*— » 

o 

£ 


or: 

3 

o 
•  « 
■*—» 

C3 

3 

cr 

p 

<4-4 

O 

£ 

P 

00 

>2 

02 


JO 

>> 

43 

60 

C 


3 

02 

£ 

P 

x: 

E- 


o 

CO 


© 

rN 


•§ 

■c 


§ 

+ 


>2 

<N 

ro 


<N 

ro 


+ 


© 


<N 

“S3 


5 

+ 

“os 


© 


<N 


H 

+ 


(N 

ro 


rs 

ro 


* 

CL 


WIT 
*  ■* 

II 

+ 

_*— » 

3 

3 

II 

>  <N 

o 

CO 

1 

1  “S3 

o 

Tf 

co 

<N 


K 

O 

3; 

+ 

© 

II 


ro 


ro 


W 

+ 


© 

«N 

^3 


0, 

"S3 


H 

+ 


rs 

ro 


0, 

ro 


.  o 
5  p 

o  x 

u-  3 
02  O 
c3  *J3 

43  g 

P  3 

3  C 

02  .O 

> _ /  «4-< 

£  •§ 

£  S3 

G  p 

P  43 

IM  ^ 
ert  rS 
cj  32 

a  p 

=  «p 

£  & 

O  P 

c  bQ 

.£ 
j-  3 
co  P 


P 

-3 

bfi  . 

•  y-  < 

^  “15 

p  ^ 
02 

P  02 
-3  c3 

4— < 

,  ^  02 

E-  4— > 

P  p 

-o  p 

•«-4  P* 
02  <4-< 

3  P 

g  P 
P  V2 

4-.  P 

O  43 
3  ~ 
_  u, 
O  p 
T3  TO 

"S3 

3 
O 
P 


<N 


£ 

O 

13 

40 


P 

3 


t-4 

o 

43 

02 

P 


*S3  ~ 

p  02 

3  4—1 

3  —f 
P  5 


3 

i  «s 


02 

P 


C 

4* 

u. 

P 

43 


W) 

#s 

‘£ 

p 


O  £ 
P  *3 

S  « 

.2  13 

33  3 

3  O 

•4— » 

O  eb 

S  M 

of 

02  4-> 

3  «■> 

3  .P 


£ 
Q-  P 
02  *£j 

3  TO 
P  cu 

E  •« 

rv  3 

o  2 

*S  ^ 

>  p 
p  p 

•o  > 

p  2 
43  ^ 
—  P 

•ft  * 


£ 


p 

4— > 
U4 

P 

*H 

o 


p 

eb 

'£ 

p 

3 


P  — 
43  W 
bD  3 

2  -B 


T3 

3 

a 

4-*t 

KA 

o 

£ 

p 


4= 

-  3 
02  O 


p 

£ 

© 

p 


02 

u 

CQ 


p 

o  H 


p 

43 


ro 


5  P  ^ 
§ 

5?  2  tS- 


o 


rs 

ro 


eN 

4»i 

ro 


+ 


© 

bL 

^3 


<N 


■S3 


5  IS 

ll 


ro 

ro 


<*2  •- 

rS 


t«4 

o 

3: 

+ 

© 


^  >2 
•—  C4 

>2  ro 


ro 


H 

S3 


51 


SYSTEM  ENERGY  &  JUDICIOUS  LIAPUNOV  FCTS. 

The  total  energy  of  the  system  (constant  in  the  absence  of  control  or  disturbances)  is: 
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Eq.  (4)  has  its  global  minimum  at  the  target  final  state: 

{M}«w  =  <0/. °).  {*<*.«),  !*«=( 0,0),  1=1,2, 3, 4 

More  generally,  error  energy  can  be  measured  from  a  time  varying  target  trajectory. 


Dissipation  of  Error  Energy  =>  Liapunov  Stable  Control  Laws 

Recall  our  candidate  Liapunov  function  is  the  error  energy  functional: 
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The  pervasive  dissipation  condition  that  U  =  -  6  is  strictly  negative ,  for  asymptotic 

stability,  is  satisfied  only  if  the  system  is  fully  controllable.  In  the  linear  case,  we  find 
that  the  anti-symmetric  in  opposition  modes,  (for  a  perfectly  symmetric  structure,  4  iden¬ 
tical  appendages)  have  zero  hub  motion  &  are  uncontrollable  by  a  hub  torque  actuator. 


Stability  Robustness:  =>  Toward  a  General  Methodology 

It  is  significant  that  if  one  includes  all  of  the  following  linear  and  nonlinear  HOT  effects: 
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Work/Energy  Principle  =>  Power  =  energy  rate  (due  to  working  forces) 
Energy  Equation  Form  Depends  upon  the  System  (or  Sub  System )  Model _ 


Weighted  Power  Decomposition  by  Substructures 
Energy  Liapunov  Functions  Stable  Control  Laws 
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Near-Minimum-Time  Maneuvers  of  DPS:  Qualitative  Issues 
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The  rigid  body  minimum-time,  bang-bang  control  can  be  generalized  to  fomally  control 
flexible  body  maneuvers,  but  we  have  found  bang-bang  control  of  flexible  body  dynamics 
usually  lacks  robustness  with  respect  to  modeling  errors .  An  approach  is  desired  which 
accommodates  the  inherent  compromise  between  minimum-time,  minimum-vibration, 
and  maximum-robustness. 


Globally  Stable  Tracking  Law 
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This  elegant  tracking  law  provides  global  Liapunov  path  stability ;  it  is  exactly  analogous 
to  the  previous  result  except  all  displacements  are  measured  from  instantaneous  point  on 
the  reference  trajectory.  The  tracking  law  is  invariant  with  respect  to  the  usual  beam 
approximations.  Potential  trouble  in  toy  land!  Minimum  time  implies  a  sense  of  urgency ! 
The  ( )r  quantities  often  cannot  be  computed  in  near  real  time. 


Stable  Tracking  With  Respect  to  an  Approximate  Trajectory 
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For  the  case  of  identical  beams  and  anti-symmetric  in  unison  deformations,  we  have  care¬ 
fully  validated  this  result  both  numerically  and  experimentally  in  our  laboratory. 


CONCLUDING  REMARKS 
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ABSTRACT 

While  excellent  progress  has  been  made  in  deriving  algorithms  that  are  efficient  for  certain  combinations  of 
system  topologies  and  concurrent  multiprocessing  hardware,  several  issues  must  be  resolved  to  incorporate 
transient  simulation  in  the  control  design  process  for  large  space  structures.  Specifically,  strategies  must  be 
developed  that  are  applicable  to  systems  with  numerous  degrees  of  freedom.  In  addition,  the  algorithms 
must  have  a  growth  potential  in  that  they  must  also  be  amenable  to  implementation  on  forthcoming 
parallel  system  architectures.  For  mechanical  system  simulation,  this  fact  implies  that 

(ii)  Algorithms  are  required  that  induce  parallelism  on  a  fine  scale,  suitable  for  the  emerging  class 
of  highly  parallel  processors. 

(iii)  Transient  simulation  methods  must  be  automatically  load  balancing  for  a  wider  collection  of 
system  topologies  and  hardware  configurations. 

This  paper  addresses  these  problems  by  employing  a  combination  range  space  /  preconditioned  conjugate 
gradient  formulation  of  multi-degree-of-freedom  dynamics.  The  method  described  herein  has  several 
advantages.  In  a  sequential  computing  environment,  the  method  has  the  features  that: 

(i)  By  employing  regular  ordering  of  the  system  connectivity  graph,  an  extremely  efficient 
preconditioner  can  be  derived  from  the  "range  space  metric",  as  opposed  to  the  system  coefficient 
matrix. 

(ii)  Because  of  the  effectiveness  of  the  preconditioner ,  preliminary  studies  indicate  that  the  method 
can  achieve  performance  rates  that  depend  linearly  upon  the  number  of  substructures,  hence  the  title 
"Order  N". 

(iii)  The  method  is  non-assembling,  i.e.,  it  does  not  require  the  assembly  of  system  mass  or 
stiffness  matrices,  and  is  consequently  amenable  to  implementation  on  workstations. 

Furthermore,  the  approach  is  promising  as  a  potential  parallel  processing  algorithm  in  that 

(iv)  The  method  exhibits  a  fine  parallel  granularity  suitable  for  a  wide  collection  of  combinations 
of  physical  system  topologies  /  computer  architectures. 

(v)  The  method  is  easily  load  balanced  among  processors,  and  does  not  rely  upon  system 
topology  to  induce  parallelism. 


(1.0)  INTRODUCTION 


There  is  no  doubt  that  an  effective  design  process  for  the  space  station  absolutely  requires 
that  high  fidelity  simulations  of  the  transient  response  to  control  inputs  be  rapidly 
attainable.  Much  research  has  been  carried  out  over  the  past  few  years  that  concentrates  on 
improving  the  performance  of  methods  for  simulating  the  dynamics  of  nonlinear, 
multibody  systems  [Gluck], [Haug], [Singh].  The  research  has  primarily  been  devoted  to 

(i)  the  derivation  of  more  efficient  formulations  of  multibody  dynamics,  and  to 

(ii)  the  derivation  of  parallel  processing  algorithms. 

Perhaps  the  most  significant  research  addressing  these  two  areas  has  been  the  introduction 
of  the  recursive,  Order  N  algorithms  in  [Hollerbach], [Feathers tone],  and  their  subsequent 
refinements  in  [Bae],  [Singh]  for  systems  of  rigid  bodies.  As  noted  in  [Singh],  these 
methods  have  the  feature  that  the  computational  cost  of  the  solution  procedure  is  linear  in 
the  number  of  degrees  of  freedom  N  of  the  system,  while  conventional  Lagrangian 
formulations  are  of  cubic  order .  The  conclusion  that  the  Lagrangian  methods  are  of  cubic 
order  derives  from  the  fact  that  a  system  generalized  mass/inertia  matrix  of  dimension  N  X 
N  must  be  factored  at  each  time  step.  Just  as  importantly,  the  computational  structure  of 
the  recursive  Order  N  algorithms  is  amenable  to  parallel  computation  for  some  system 
topologies.  If  the  system  to  be  modelled  has  many  independent  branches  in  its  system 
connectivity  graph,  the  computational  work  required  by  the  algorithm  can  be  distributed 
among  processors  by  assigning  branches  to  independent  processors.  As  an  example, 
figure  (1.1)  illustrates  the  connectivity  graph  for  an  all  terrain  vehicle  modelled  in  [Bae]. 
The  processor  computational  load  distribution  employed  in  the  paper  is  likewise  illustrated. 
Because  of  the  system  connectivity  and  specific  hardware  architecture,  excellent 
performance  improvements  and  processor  utilization  are  achieved  in  [Bae]. 

Due  to  these  successes  for  rigid  body  simulations,  it  is  well-known  that  many  research 
institutions  are  presently  investigating  adaptations  of  the  original  recursive  method  to  model 
systems  comprised  of  flexible  bodies.  No  doubt,  the  result  will  be  highly  efficient 
algorithms  that  perform  well.  Still,  three  key  goals  must  be  resolved  before  a  general 
parallel  processing  algorithm  can  be  obtained. 

(i)  Algorithms  are  required  that  induce  parallelism  on  a  finer  scale,  suitable  for  the 
emerging  class  of  highly  parallel  processors. 

(ii)  Concurrent  transient  simulation  methods  must  be  automatically  load  balancing 
for  a  wider  collection  of  combinations  of  mechanical  systems  and  concurrent 
multiprocessing  hardware. 

(iii)  The  transient  simulation  method  should  also  be  amenable  to  vector  processing 
implementation  on  each  independent  concurrent  multiprocessor. 

Based  upon  preliminary  investigation,  these  goals  should  be  very  challenging  if  the 
algorithm  is  based  upon  an  recursive  Order  N  formulation. 
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FIGURE  (1.1)  -  PROCESSOR  ASSIGNMENT 


An  innovative  strategy  based  upon  these  goals  is  derived  in  this  paper.  In  part,  its 
foundation  can  be  traced  to  element-by-element  methods  already  in  use  in  finite  element 
solution  procedures  [Hughes].  As  regards  sequential  computing  environments: 

(i)  The  combination  range  space  formulation  /  PCG  solution  is  an  extremely 
efficient  sequential  algorithm  for  a  class  of  problems  described  in  the  paper.  The 
efficiency  is  primarily  due  to  the  selection  of  a  Block  Jacobi  preconditioner  that  is 
rapidly  convergent. 

(i)  The  method  is  non-assembling,  i.e.  ,  it  does  not  require  a  large  amount  of  in- 
core  storage,  and  consequently  is  also  attractive  as  a  candidate  for  implementation 
on  workstations. 

(iii)  Preliminary  studies  indicate  that  due  to  the  rapid  convergence  achieved  by 
using  the  selected  preconditioner,  the  method  can  achieve  performance  rates  that 
depend  linearly  upon  the  number  of  substructures. 

Moreover,  the  method  should  be  readily  implemented  on  parallel  processors: 

(iii)  A  vast  literature  exists  on  the  amenability  of  the  PCG  solution  procedure  to 
both  concurrent  and  vector  processing. 

(iv)  The  method  is  relatively  easily  load  balanced  among  processors,  and  does  not 
rely  upon  system  topology  to  induce  parallelism. 

This  paper  focuses  on  the  fundamental  dynamical  formulation  using  a  combination  range 
space  /  PCG  solution,  and  its  performance  on  sequential  computing  machines.  Although 
the  potential  applicatior  of  the  method  on  parallel  architectures  is  outlined,  the  details  of  a 
concurrent  implementation  are  presented  in  a  forthcoming  paper. 
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(2.1)  RANGE  SPACE  /  PRECONDITIONED  CG  EQUATIONS 


The  range  space  formulation  of  dynamics  has  been  derived  in  the  aerospace  and  mechanism 
dynamics  research  literature  in  [Placek],[Agrawal],[Kurdila].  Its  theoretical  foundation  can 
be  traced  to  the  range  space  formulation  of  constrained  quadratic  optimization  [Gill].  Still, 
despite  the  fact  that  it  is  often  less  computationally  expensive  than  the  nullspace  methods, 
the  nullspace  method  seems  to  have  received  more  attention  in  the  literature  [Singh, 
Wehage,  Kim,  Huston,  Kurdila...].  If  the  dynamics  of  a  nonlinear,  multibody  system  are 
governed  by  the  collection  of  differential-algebraic  equations 

[M(q)q=  f(q,q,t)  +  [C(q)]TX 

subject  to  constraints  in  linear,  non-holonomic  form 


[C[q]]q=  0 

the  range  space  solution  of  these  equations  are  given  by  explicitly  solving  for  the 
multipliers 

*=-([C(<J)][M(<J)r'[C(<j)f  )'’{[C(«j)I*»(<J)rV  «7,<j  ,  0  -  e(q,  q,  ()} 

and  substituting  to  achieve  a  govering  system  of  ordinary  differential  equations. 

q  =  [M(q)]\f  ( q .  q .  0  -[£(<?)  ]r 

{([C(iJ)][M(<J)]‘'[C«J)]r)''{[C(q)IM«))]'V  (q,q  ,  t)  -  e(q,  q,  »)}} 

In  the  above  equations,  the  constraints  have  been  differentiated  twice  to  yield 


[C(q)]q  =  -  -~j~j([C{q)])q  =  e{q,  q,  t ) 

Any  standard  explict-predictor  /  implicit-corrector,  or  Runge-Kutta  integration  scheme  can 
be  applied  to  these  equations  provided  that  the  condition  number 


KJA) 


of  the  constraint  metric 

[C(q)lM(q)]'\c(q)]T 

does  not  become  too  large.  The  restriction  that  the  condition  number  above  remains  small 
precludes  the  possibility  of  redundant  constraints  (for  example,  as  associated  with 
singularities  arising  from  closed  loops)  and  remains  an  underlying  assumption  throughout 
the  rest  of  the  paper. 
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One  advantage  of  the  range  space  equations  for  systems  having  many  independent 
structures  to  be  assembled  is  that  the  system  coefficient  matrix  is  block  diagonal  and, 
consequently,  the  factorization  and  back-substitution  required  to  form  the  product  of  the 
inverse  of  the  mass  matrix  and  a  given  vector  is  relatively  inexpensive  to  calculate.  It 
requires  that  one  calculate  the  factorization  of  the  individual  substructure  mass  matrices 
alone.  In  fact,  one  need  not  even  assemble  the  system  mass  matrix,  and  the  factorizations 
can  occur  in  parallel.  Unfortunately,  if  one  subdivides  the  overall  system  into  finer 
collections  of  substructures  (to  facilitate  the  factorization  of  the  system  coefficient  matrix), 
numerous  constraints  are  introduced  into  the  model.  As  a  consequence,  one  has  the 
tradeoff  shown  in  figure  (2.1.1). 

The  approach  taken  in  this  paper  is  to  finely  subdivide  the  system  to  be  modelled,  and  thus 
accrue  the  benefits  of  having  a  system  coefficient  matrix  with  smaller  block  diagonals,  but 
also  employ  a  solution  procedure  that  ameliorates  the  cost  associated  with  the  increasing 
dimensionality  of  the  constraint  metric.  Specifically,  the  calculation  of  the  Lagrange 
multipliers  in 

A=-([C(<j)][/W(<j)r'[C«,)]r)‘'{[C(<7)l[M(q)fV  (q,q  ,  t )  -  e«j,  q,  I)} 

is  carried  out  using  the  preconditioned  conjugate  gradient  procedure. 

(2.2)  THE  PRECONDITIONED  CONJUGATE  GRADIENT  SOLUTION 

The  preconditioned  conjugate  gradient  procedure  is  an  "accelerated"  variant  of  the  classical 
conjugate  gradient  procedure.  If  it  is  required  to  solve  the  linear  system  of  equations 

Ax  =  b 

the  procedure  can  be  summarized  from  [Golub] 

*o=0 

ro=b 

For  k  =  1, ...  n 

then 

x  =  x 

*  Ak- 1 

else 

Solve  Ozk_=rk  i 

=  r,-,  -  a,APk 
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FIGURE  (2.1.1)  :  COMPUTATATIONAL  TRADEOFF  DUE  TO 
SUB  DOMAIN  DECOMPOSITION 


Careful  inspection  of  the  algorithm  shows  that  the  most  computationally  expensive  tasks  in 
the  procedure  are  the 


(i)  calculation  of  the  product  of  the  coefficient  matrix  A  and  a  given  residual 
vector, 


(ii)  and  the  solution  of  a  linear  system  of  equations  requiring  the  factorization  of  the 
preconditioner  Q  . 

The  rate  of  convergence  of  the  preconditioned  conjugate  gradient  algorithm  is  accelerated 
by  employing  a  user-defined  "preconditioning  matrix."  This  matrix  must  have  two 
properties  to  be  an  effective  preconditioner: 

(i)  It  must  be  relatively  easy  to  factor. 

(ii)  It  must  be  an  approximate  inverse  to  the  constraint  metric  in  a  sense  to  be  made 
precise  below. 


The  reason  for  employing  the  preconditioned  conjugate  gradient  solution  method  is  that  the 
convergence  rate  of  the  conjugate  gradient  algorithm  (that  is,  with  Q  =  I)  is  governed  by 


< 


~1  -V^(A) 

.1  +  Vic  (A) 


Thus,  the  rate  of  convergence  of  the  algorithm  improves  as  the  condition  number 

k(A) 


decreases.  The  reference  [Golub]  has  shown  that  the  convergence  of  the  preconditioned 
conjugate  gradient  method  is  governed  by  the  same  expression,  but  with  A  replaced  with 

a=  q~2aqj 

Clearly,  if  the  preconditioner  is  identical  to  the  coefficient  matrix,  then  the  condition 

number  of  ^  is  minimized.  Hence,  the  preconditioner  is  sought  such  that  its  inverse 
approximates  the  inverse  of  the  coefficient  matrix.  Many  methods  exist  for  the  calculation 
of  preconditioners  [Golub].  It  should  be  noted  that  while  the  motivation  for  the  use  of 
many  of  these  preconditioners  is  mathematically  sound,  the  final  choice  invariably  involves 
some  heuristic. 
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(2.3)  THE  CHOICE  OF  THE  PRECONDITIONER 


The  choice  of  the  preconditioner  employed  in  this  paper  is  based  upon  the  following 
assumptions  regarding  the  structural/mechanical  system  to  be  modelled: 

(i)  The  system  closely  resembles  a  series  of  chains  of  bodies 

(ii)  The  number  of  interface  degrees  of  freedom  is  small  relative  to  the  number  of 
interior  degrees  of  freedom  for  a  substructure. 

(iii)  The  system  does  not  contain  any  closed  chains. 

To  a  large  extent,  these  assumptions  have  been  driven  by  the  physical  structure  of  the  space 
station  in  its  assembly  complete  configuration. 

The  preconditioner  for  the  system  constraint  metric  is  based  upon  the  topology  of  a  chain  of 
substructures  as  shown  in  figure  (2.3.1).  If 

d  xN 

denotes  the  constraint  matrix  connecting  two  bodies  at  the  ith  interface,  the  system 
constraint  matrix  has  the  form 


[C«7)]  = 


[C2] 


Ff'" 


The  system  constraint  metric  can  then  be  written 


\c)[M]~\c)T  -  [c}[M]-\ck]T 

[c2][M]-\cf 

_[ck][M]\cf  ...  [ck][M)\ck]\ 

Based  upon  the  structure  of  the  constraint  metric  above,  the  preconditioner  is  selected  to  be 
the  block  diagonal  matrix 
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FIGURE  (2.3.1)  :  MULTIBODY  CHAIN 


[c2][M]'\c2]t 

[ck][M]~\ck]T 

Although  the  off-diagonal  blocks 

[Qmwm'iCiiq)]7  -  o 

(for  i  not  equal  to  j)  are  not  generally  identically  equal  to  zero,  this  choice  of  preconditioner 
is  shown  to  be  extremely  efficient  for  the  class  of  problems  described  in  the  next  section. 
Furthermore,  this  preconditioner  satisfies  the  two  essential  criteria  of  good  preconditioners: 

(i)  It  is  block  diagonal,  with  small  diagonal  blocks,  and  is  relatively  easy  to  factor. 

(ii)  It  has  an  inverse  that  provides  a  good  approximation  to  the  inverse  of  the  full 
system  coefficient  matrix. 

This  latter  conclusion  results  from  the  well-known  fact  [Wittenburg]  that  the  directed  graph 
representing  the  connectivity  of  an  open  loop  system  can  be  regularly  ordered.  The  regular 
ordering  results  in  a  system  constraint  metric  that  has  a  reduced  bandwidth.  That  is,  many 
of  the  off-diagonal  blocks 

[C«7)][M(c7)]'1[C/((7)]r  =  0 

are  identically  zero  for  i»j.  The  choice  of  preconditioner  shown  above  is  often  denoted 
the  Block  Jacobi  preconditioner  and  is  known  to  be  highly  effective  for  classes  of  systems 
of  equations  arising  from  elliptic  partial  differential  equations  [Reid], 
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(3.1)  SPACE  MAST  SIMULATIONS 


To  establish  the  efficiency  and  performance  of  the  combined  range  space  /  PCG  algorithm, 
several  transient  simulations  have  been  carreid  out.  The  first  such  simulation  has  been 
designed  to  answer  two  key  questions  regarding  the  feasibility  of  the  approach  for 
multibody  structures: 

(i)  How  efficient  is  the  selected  preconditioner  and  how  does  the  preconditioner 
improve  the  effectiveness  of  the  conjugate  gradient  iteration  processes  as  a  whole? 

(ii)  How  efficiently  can  the  range  space  /  PCG  solution  for  the  transient  response 
be  carried  out?  In  particular,  what  is  the  computational  cost  of  solving  the 
"constraint  metric"  factorization  at  each  time  step? 

These  two  questions,  that  is,  preconditioner  efficiency  and  constraint  metric  factorization 
were  judged  to  the  crucial  computational  "bottlenecks"  for  the  algorithm  as  a  whole. 

The  two  substructures  selected  for  evaluation  of  the  algorithm  in  studying  its  feasibility  are 
shown  in  figures  (3.1.1)  and  (3.2.2).  The  first  structure  is  a  63  degree  of  freedom  Z-truss 
substructure  comprised  of  rod  elements  generated  by  MSC  PAL.  The  second  assembly  is  a 
54  degree  of  freedom  X-ffame  substructure  comprised  of  three  dimensional  beam  elements 
generated  by  the  same  program.  Figures  (3.1.3)  and  (3.1.4)  illustrate  that  a  remarkable 
convergence  rate  is  achieved  using  the  preconditioner  derived  from  the  regularly  ordered 
graph  of  the  constraint  metric.  In  figure  (3.1.3)  the  number  of  iterations  required  by  the 
preconditioned  conjugate  gradient  algorithm  to  converge  to  a  tolerance  of  l.e-14  of  the  true 
solution  is  plotted  versus  the  number  of  degrees  of  freedom  on  the  horizontal  axis.  Thus, 
the  number  of  flexible  degrees  of  freedom  on  the  horizontal  axis  varies  from  126  to  504  as 
substructures  are  contilevered  end-to-end.  As  clearly  illustrated  in  the  diagram,  the  number 
of  preconditioned  conjugate  gradient  iterations  remains  constant,  and  equal  to  3, 
independent  of  the  number  of  total  degrees  of  freedom.  The  other  line  plotted  on  the  graph 
is  the  analytical  upper  limit  to  the  number  of  iterations  required  in  the  conjugate  gradient 
method.  Completely  analogous  results  are  depicted  in  figure  (3.1.4).  In  this  case  the 
number  of  total  flexible  degrees  of  freedom  in  the  structure  varies  from  108  to  324,  and  the 
number  of  iterations  of  the  preconditioned  conjugate  gradient  algorithm  required  to  achieve 
convergence  to  a  tolerance  of  l.e-14  is  plotted  on  the  vertical  axis.  Again,  independent  of 
the  number  of  degrees  of  freedom,  the  number  of  iterations  remains  constant  and  equal  to 
4.  Consequently,  one  can  conclude  that  the  algorithm  for  generating  the  preconditioner  is 
indeed  extremely  efficient,  for  the  test  problems  simulated. 

Figures  (3.1.5)  and  (3.1.6)  depict  the  total  time  required  per  time  step  to  solve  the 
"inversion"  of  the  system  constraint  metric.  Figure  (3.1.6)  depicts  the  corresponding 
results  for  the  Z-truss,  while  figure  (3.1.6)  depicts  the  results  for  the  X-frame.  In  both 
cases  it  is  clear  that  the  simulation  grows  linearly  as  a  function  of  the  total  number  of 
flexible  degrees  of  freedom.  To  the  authors'  knowledge,  no  such  result  for  flexible  bodies 
has  been  presented  to  date. 

In  later  sections,  while  considering  the  amenability  of  the  overall  solution  procedure  to 
parallel  processing,  it  is  crucial  to  consider  where  computation  time  is  spent  during  a  typical 
time  step.  These  results  are  shown  in  figure  (3.1.7).  This  figure  shows  that  two  steps 
dominate  the  time  spent  solving  the  range  space  /  PCG  solution  procedure: 

(i)  the  time  spent  forming  the  product  of  the  inverse  of  the  system  mass  matrix  and 
a  given  vector,  and 
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(ii)  the  time  spent  applying  the  inverse  of  the  preconditioner  to  a  given  vector. 

As  will  be  discussed  in  more  detail  in  a  later  secdon  on  applications  in  parallel  processing, 
these  two  steps  are  easily  parallelized  because  of  their  block  diagonal  structure,  and  should 
yield  excellent  improvements  in  performance  on  concurrent  multiprocessors. 


Z-TRUSS  SUBSTRUCTURE 
MSC/PAL  ROD  ELEMENT 
63  DOF  PER  SUBSTRUCTURE 


Z-TRUSS  SUBSTRUCTURE 
NASTRAN  ELEMENT  TYPE  2 
54  DOF  PER  SUBSTRUCTURE 


FIGURE  (3.1.2)  :54  DOF  X-FRAME  SUBSTRUCTURE 
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Substructure  2 


Substructure  1 


NUMBER  OF  ITERATIONS 


ITERATIONS  AS  A  FUNCTION  OF  SUBSTRUCTURES 


FIGURE  (3.1.3) : ITERATIONS  FOR  Z-TRUSS 


ITERATIONS 


ITERATIONS  AS  A  FUNCTION  OF  SUBSTRUCTURES 


NUMBER  OF  BAYS  (DOF/54  :  54  DOF  PER  SUBSTRUCTURE) 


FIGURE  (3.1.4) : ITERATIONS  FOR  X-TRUSS 


79 


DOF  (! 


FIGURE  (3.1.5)  :SYSTEM  CONSTRAIN 


Range  Space  Metric  Factorization  -  Breakdown 


FIGURE  (3.1.7)  .BREAKDOWN  OF  TIMING  RESULTS  FOR  63  DOF  Z-TRUSS 


x(74)  (Inches) 
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PCG  Iterations 


Time  Integration  (1008  DOF,  0(h**4)  RK) 


FIGURE  (3.1.9) :  PCG  ITERATION  TIME  HISTORY  /  63  DOF  Z-TRUSS  /  1008  DOF  MODEL 
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(3.2)  SPACE  STATION  SIMULATIONS 


While  the  simulations  of  the  truss  structures  in  the  last  section  answer  important  questions 
regarding  the  efficiency  and  performance  of  the  PCG  algorithm,  it  remains  to  be  shown  that 
the  method  is  robust  in  problems  that  do  not  enjoy  such  structure  in  form.  By  analogy,  one 
can  liken  the  addition  of  identical  substructures  with  identical  interfaces  to  the  solution  of 
Poisson's  equation  on  a  rectangular  grid:  many  methods  exist  that  will  have  "remarkable" 
convergence  rates  for  such  highly  structured  problems.  Unfortunately,  the  performance 
degrades  rapidly  as  more  general  problems  are  considered. 

In  this  section,  the  simulation  of  Space  Station  Freedom  in  the  assembly  complete 
configuration  is  considered.  Qualitatively,  the  simulations  described  in  this  section  differ 
from  those  in  the  last  section  in  three  important  respects: 

(i)  the  number  of  degrees  of  freedom  per  substructure  varies, 

(ii)  the  constituent  substructures  matrices  are  not  identical  in  general  (although 

some  occur  in  symmetric  pairs  about  the  core  body), 

(iii)  the  number  of  degrees  of  freedom  between  substructure  interfaces  varies. 

The  model  employed  is  a  simplified  version  of  the  full  finite  element  model  for  the 
assembly  complete  space  station  shown  in  figure  (3.2.1).  The  space  station  is  subdivided 
into  13  individual  substructures  having  from  90  to  150  degrees  of  freedom  each.  The 
number  of  constraints  per  interface  varies  from  24  to  36.  The  first  simulation  considers 
only  the  central  bodies  5  through  9.  Blowups  of  these  five  bodies  are  shown  in  figures 
(3.2.2)  through  (3.2.6),  and  their  position  in  the  assembly  complete  station  is  shown  in 
figure  (3.2.7).  A  concise  summary  of  the  results  of  several  simulations  is  given  in  figures 
(3.2.8)  and  (3.2.9).  Figure  (3.2.9)  plots  the  total  time  per  integration  time  step  versus  the 
number  of  degrees  of  freedom  in  the  model.  The  model  is  assembled  from  left  to  right,  so 
that  data  points  are  plotted  for  systems  comprised  of  2,  3,  4,  and  5  substructures.  The 
largest  model  considered  is  comprised  of  nearly  500  degrees  of  freedom.  Because  the 
number  of  degrees  of  freedom  per  body  and  number  of  constraints  per  interface  varies  with 
the  addition  of  each  set  of  bays,  the  computational  cost  plots  are  not  exactly  linear  as  in  the 
space  boom  simulations.  Still,  a  linear  curve  fit  describes  the  timing  data  well  for  all  the 
simulations  considered.  As  in  the  case  of  the  space  boom  simulations,  a  very  large  fraction 
of  the  overall  simulation  time  is  spent 

(i)  calculating  a  product  of  the  inverse  of  the  mass  matrix  and  a  given  vector, 

(ii)  calculating  the  product  of  the  inverse  of  the  preconditioner  and  a  given  vector. 

As  noted  earlier,  both  of  these  matrices  are  block  diagonal  and  trivially  parallelizable. 
Hence,  the  potential  for  implementing  the  method  in  a  parallel  environment  is  very 
promising.  Figure  (3.2.10)  depicts  the  number  of  PCG  iterations  versus  the  integration 
time  history.  As  in  the  previous  example,  the  selected  preconditioner  is  extremely 
effective,  and  requires  4  iterations  for  convergence,  essentially  independent  of  the  number 
of  degrees  of  freedom  in  the  model. 
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FIGURE  (3.2.2) :  SUBSTRUCTURE  5 


FIGURE  (3.2.3) :  SUBSTRUCTURE  6 
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FIGURE  (3.2.4) :  SUBSTRUCTURE  7 
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FIGURE  (3.2.5) :  SUBSTRUCTURE  8 


FIGURE  (3.2.6) :  SUBSTRUCTURE  * 


FIGURE  (3.2.7) :  SUBSTRUCTURE  LOCATIONS  IN  ASSEMBLY  COMPLETE  CONFIGURATIONS 
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Space  Station  Model  -  5  Substructure  Breakdown 


FIGURE  (3.2.8) :  SPACE  STATION  MODEL  -  5  SUBSTRUCTURE  TIMING  BREAKDOWN 
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PCG  Iterations 


Space  Station  /  5  Body  Model  /  PCG  Iterations 
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Tims  (x  .0001  sec) 


FIGURE  (3.2.9) :  PRECONDITIONED  CONJUGATE  GRADIENT  ITERATION  TIME  HISTORY 
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ABSTRACT 

While  earlier  papers  have  studied  the  convergence  properties  of  Lyapunov  stable  penalty  methods 
(LSPM)  as  applied  to  spectral  approximations,  the  current  paper  investigates  the  existence  and  con¬ 
vergence  of  penalty  approximations  when  applied  to  transient  analysis.  This  paper  makes  use  of 
standard  techniques  in  the  analysis  of  linear,  hyperbolic  partial  differential  equations  to  show  that 
a  sequence  of  solutions  generated  by  the  Lyapunov  stable  penalty  equations  approaches  the  solu¬ 
tion  of  the  differential-algebraic  equations  (DAE’s)  governing  the  dynamics  of  multibody  prob¬ 
lems  arising  in  linear  vibrations.  Specifically,  the  analysis  relies  upon 

(1)  the  existence  of  Lyapunov  functions  for  the  class  of  problems  considered, 

(2)  standard  operator  norm  approximations  of  orthogonal  projections  onto  the  range  of  the 
transposed  constraint  matrix, 

(3)  the  application  of  Gron wall’s  inequality  to  show  that  the  sequence  of  approximate  so¬ 
lutions  remains  bounded  for  compact  intervals  in  time. 

The  analysis  is  quite  general  in  that  no  assumption  is  made  that  the  system  be  natural  or  completely 
integrable.  The  result  of  the  analysis  is  the  derivation  of  an  explicit  variational  relationship  between 
the  norm  of  the  constraint  violation  time  history  and  the  error  between  the  solutions  of  the  true  and 
approximate  penalty  equations.  For  linear,  undamped  multi-degree-of-freedom  equations  this  re¬ 
lationship  takes  the  form 


{ |j  *  (*  -  *£)  || }  *  0||  /*^e||  *  a,  |  *e || 


In  the  case  of  damped  multi-degree-of-freedom  equations,  the  variational  relationship  is 


(1.0)  INTRODUCTION 


The  difficulty  that  may  arise  in  numerically  integrating  systems  of  differential-algebraic  equations 
is  well-documented  [Gear],[Wehage].  Numerical  procedures  intentionally  designed  for  as  general 
classes  of  DAE’s  have  been  the  topic  of  much  research  over  the  past  two  years  by  numerical  ana¬ 
lysts  [Gear],[Burrage],[Alexander]. 

In  computational  mechanics,  on  the  other  hand,  methods  have  arisen  that  utilize  the  fact  that  the 
particular  system  of  differential-algebraic  equations  to  be  solved  have  been  derived  from  Lagran- 
giar.  or  Hamiltonian  formulations  of  dynamics.  Examples  of  this  type  of  approach  are  given  in 
[Bayo  1,2,3],  [Park],  [Kurdila].  Essentially,  all  of  these  methods  approximate  the  system  of  gov¬ 
erning  differential-algebraic  equations  by  an  altogether  different  dynamical  system;  one  that  has 
been  obtained  via  penalty  perturbation  of  the  Lagrangian  or  Hamiltonian  for  the  system. 

While  all  of  thes  papers  present  considerable  empirical  evidence  that  the  penalty  methods  are  stable 
and  convergent,  little  analysis  has  been  conducted  to  establish  this  fact.  Because  the  governing 
equations  are  neither  linear,  nor  coercive  in  general,  standard  results  as  in  [Oden  1]  or  [Oden  2]  are 
not  directly  applicable.  Simple,  but  effective  error  estimates  for  spectral  approximations  using  the 
penalty  method  have  been  presented  in  [Kurdila].  Nonlinear  stability  and  convergence  criteria  are 
considered  in  [Kurdila  1],  [Kurdila  2],  but  rely  upon  the  restrictions  that  the  governing  system  be 
natural,  and  in  some  cases  conservative.  In  the  latter  case,  much  of  the  arguments  presented  are 
based  upon  the  underlying  Hamiltonian  structure  of  the  systems  considered. 

The  purpose  of  this  paper  is  to  investigate  convergence  criteria  for  linear  multi-degree-of -freedom 
systems  arising  in  linear  vibrations.  The  analysis  that  follows  is  quite  general  in  that  it  does  not 
require  that  the  system  be  conservative.  No  restriction  on  the  form  of  the  forcing  term  is  enforced 
other  than  it  is  Lj-integrable  on  bounded  intervals  of  time.  The  paper  concludes  by  deriving  vari¬ 
ational  statements  that  bound  the  error  in  approximation  by  the  norm  of  the  constraint  violation  ob¬ 
tained  in  the  approximate  solutions.  These  variational  statements  are  of  great  practical  importance: 
they  imply  that  by  monitoring  the  constraint  violation  one  can  be  assured  that  the  solution  is  accu 
rate.  One  should  note  that  this  result  is  not  true  in  general  for  nonlinear  systems!  For  the  nonlinear 
case,  bifurcations  as  described  in  l Kurdila]  can  occur  in  which  the  constraint  violation  remains 
small,  but  the  penalty  solution  diverges  from  the  true  solution. 
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(2.0)  LINEAR,  UNDAMPED  EQUATIONS  OF  MOTION 


(2.1)  Original  Equations 

As  a  starting  point  for  the  analysis,  we  consider  the  system  of  undamped,  second  order  ordinary 
differential  equations  subject  to  holonomic  constraints 

rd<Dir,  - 
Mx  +  Kx  =  =r -  X  +  F 

Ld*  J 

where 


and  the  dimensions  of  the  constituent  matrices  and  vectors  are 


xeRN 
<De  Rd 
XeRD 


e  R 


OxN 


Fe  RN 
Me  RNxN 
K  e  rn*n 


As  is  usually  encountered  in  linear  vibration  equations,  M  is  symmetric,  positive  definite,  and  may 
be  considered  to  have  been  generated  by  consistent  finite  element  formulation  (as  opposed  to  a 
lumped  formulation).  The  stiffness  matrix  K  is  assumed  to  be  symmetric  positive  semi-definite, 
and  the  constraint  matrix  is  constant.  To  simplifiy  the  analysis,  we  introduce  the  following  change 
of  variables 

.-1/2  .-1/2 
K  =  M  KM 
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.  1/2 

X  -  M  x 


.-1/2. 

F  =  M  F 

With  these  definitions,  the  governing  equations  take  the  simple  form 


X  +  KX 


T 

X  +  F 


No  generality  is  lost  in  assuming  the  equations  have  this  form,  while  considerable  simplification  is 
achieved  in  the  derivations  that  follow. 

By  differentiating  the  constraints,  and  defining  the  orthogonal  projection  P  onto  the  range  of  the 
transposed  constraint  matrix 
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the  governing  equations  can  be  expressed  in  constraint-reactionless  form 

X  +  (/- P)  KX  =  (I-P)F 
X  +  QKX  =  QF 


In  the  above  equations. 


Q  =  I-P 


is  the  orthogonal  projection  onto  the  space  of  admissable  configurations. 

(2.2)  Undamped  Penalty  Equations 

While  several  forms  of  the  penalty  equations  have  appeared  in  the  literature  [Arnold],  [Bayo  1,2,3], 
[Kurdila  1,2]  and  [Park  1],  the  form  chosen  here  can  be  derived  as  in  [Bayo]  or  [Kurdila]  from  the 
penalized  kinetic  and  potential  energies 


2  2  e 


2  2  e 


The  penalty  form  of  Lagrange’s  equations 


Lt  =  Tt~V£ 


dL.  \  dL. 


d(d±A. 

dt\dX  ) 


then  result  in 


ir3<I>ir  .. 

Xt  +  KX£  =  {O  +  aO}  +F 
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Xt  +  KXc 


1 

"dO- 

T  | 

"30] 

y  a  r9(I>l 

rrdO- 

£ 

_axj 

I 

Lax, 

Xe  el_3X. 

1  Lax. 

',  +  e 


■aon5 

Yao- 

a  raon1 

"ao- 

.3X. 

.ax. 

}^+t^+j[s 

.ax. 

}X,_  =  F 


It  is  easy  to  verify  that  f  Albert] 


-ao- 

rao- 

.ax. 

Lax. 

is  invertible  for  every  e>0,  and  that 


rao/ 

raoi 

aoi 

T  | 

-ao- 

LaxJ 

Lax, 

}  =  e  {e  + 

.ax. 

_3x_ 

-l 


This  identity  allows  one  to  rewrite  the  governing  equations 


x+f  eL5*J  ISxJ1  '  e 


ao 

axj 


rdO 


axJ 


1 


ao 


*  r90 

{/+  ilaxj  Lax 


in  the  equivalent  form 


v  r  ^aOi 

x  +  e{e  +  Uxj 


do 

LaxJ 


-I 


do 


>  <X+“[H 


H]>x<= 


rao 


£{£+  5v 


LdXj 


-ao 

.3X_ 


I 


By  carrying  out  the  multiplications 


the  equations  can  be  written  in  symbolic  form  as 

X£  +  Q(c)KX£  +  aP(e)  =  Q(e)F 

where  we  have  introduced  the  definitions 


One  should  note  that 

I-P(e)  =  Q{z) 

and  that 

P->P(e) 

Q->Q(e) 

in  the  bounded  linear  operator  topology. 
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(2.3)  Boundedness  of  Penalized  Solutions  For  Undamped  Case 

Before  establishing  the  convergence  of  the  penalty  approximations,  it  is  first  necessary  to  show  that 
the  approximations  remain  bounded 


Fs(7')||*C1 

\\x£(T)  II  <c2 


for  any  arbitrary,  fixed  final  time  T.  To  this  end  one  can  take  the  inner  product  of  the  penalized 
governing  equations  with  the  derivative  of  the  penalized  solution  to  obtain 


1  r3cb"irr30i  ..  1  rd<PY  ^cKPi  .. 

«'+EL3xj  Ls*],x-*>-  <**■> 


1 


This  expression  can  be  re-written  in  a  time  rate  of  change  of  energy  form: 


4{1| 
dr  v 


+  (KXt,Xe)}  + 
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A  single  integration  in  time  yields 


\\\x£(t)  ||2+  (KXc(t),Xe(t)) 


(  1  .7 

r 

rao- 

+  ! 

[jiX'  1 

Lax. 

Lax. 

a 


)*E(0  +2ixt«)  ( 


-aon 

ao- 

.ax! 

_axJ 

)X£(0 

y 

=  \  ||XE  (0)  ||2  +  (KXt  (0) ,  Xe (0) > 


i  .r  rd®iTrdo 

S*e<0)<L3xJ 


n  a  t  ra<Jnrra<l>'i 

])*f(0)+^<0)([^J  L3x])X«(m 

Jo<F(t),  Xt(x))dx 


dX 

+ 


But  since  all  penalized  equations  are  required  to  satisfy  the  constraints, 


a<u 


M  iXe(°)  =  |_5xjXe(0>  =  0 


the  energy  expression  becomes 


i||XE(0  \\2+(KXz(t),Xz{t)) 


+  (s^  <[jx]  [»]>*■“> +  s^w< 

=  ij|XE  (0)  ||2+<XX£(0),X£(0)> 


-ao- 

-ao- 

_ax_ 

_ax. 

)XAt) 


+  J'0<F(t),Xe(t)  )dx 
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In  particular,  this  implies  that  we  can  write  the  inequality 


^lA'e(r)  ||2+<KXE(t),Xe(0> 
<^(0)||2+(tfX£(0),X£(0)> 
+  Joll^(t)  li ||Xe(x)  | dx 


But  by  using  the  (trivial)  form  of  Minkowski’s  inequality  [Oden] 

2\ab\  <a2  +  b2 


one  can  write 


~;XE(r)jl2  +  <KXE(r),XE(r)> 

<5:iX£(0)j!2  +  ~(0E(0)>XE(0)) 

+  F(T)  :i2dT+^J‘l|XE(T)  fdX 

~l\Xz(t)  ||2+l(XXE(t),XE(0) 

<  {  i  "\Xt  (0)  |j2  +  l-  (KXt  (0) ,  X£  (0) )  +  I  Jill  F  (x)  ||2dx} 
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Finally,  one  can  write  the  energy  expression  in  a  form  appropriate  for  the  application  of  Gronwall’s 
inequality  [Wloka] 


i||*£(r)  ||2+l<KX£(r),X£(r)> 

<  {it|Xe(0)|2+l  (KXt  (0) , XE  (0)  >  +  \  JiH F  (x)  || 2dx} 
+  Ij!>(||Xe(x)  ||2+(/fX£(T),X£(x)»cfx 


~,Xt  it)  2  +  \  (KXt  (t),X£(t))<g  (t)  +  \  j‘Q  (||X£  (x)  ||2  +  1  (KXt  (x) ,  XE  (x) ))  dx 


Application  of  Gronwall’s  inequality  to  the  above  expression  implies  the  desired  result  that  the  pe¬ 
nalized  solution  remains  bounded 


I Xe(0  |  <K, 
(KXE(t),XE(t))<  k2 


for  a  fixed  time  t,  0<t<T. 

(2.4)  Variational  Dependence  of  Error  on  Constraint  Violation 

With  the  above  result  implying  that  the  solution  of  the  penalty  equations  remains  bounded,  one  can 
now  derive  a  relationship  between  the  error  in  the  penalty  approximation 

\\x-x£\\ 
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and  the  norm  of  the  constraint  violation 


1 

racn 

-l 

||  2  = 

M. 

_dX_ 

) 

Subtracting  the  original,  exact  equations  and  the  penalty  form  of  the  equations,  one  can  obtain 

X-XE  +  QKX-Q(e)  kxe 
=  ( Q-Q(e))F  +  aP(e)Xt 


By  the  addition  and  subtraction  of  identical  terms,  and  by  employing  the  triangle  inequality,  one 
can  write  that 


I!  X  -  X£|  +  ||  Q  -  Q  (£)  ||  ||  KXt\\  +  II  Q\\  ||  K  (X  -  XE)  || 
^  ;|<2-<2(e)  IIHFil  -Hoc|jPATe||  +a||P  (e) -P\\  ||Xe 


Using  the  uniform  convergence  of  P(e)  to  P,  the  fact  that  the  solution  of  the  penalty  equations  re¬ 
mains  bounded,  the  final  variational  error  inequality  is  achieved. 
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(3.0)  LINEAR,  DAMPED  MDOF  EQUATIONS 
(3.1)  Original,  Exact  Equations 

The  derivation  of  a  corresponding  variational  statement  for  the  case  of  damped,  MDOF  systems 
follows  in  much  the  same  manner  as  the  strategy  empolyed  in  the  undamped  case.  In  this  section, 
the  governing  system  of  differential-algebraic  equations  are 


Mx  +  Cx  +  Kx 


=  0 


where,  as  in  the  undamped  case, 


xeRN 
Oe  R° 

Xe  Rd 

\  P 1  s  Rd*n 


FeRN 
Me  RN*N 
Ke  RNyN 
C  6 


The  same  assumptions  regarding  the  properties  of  the  mass,  stiffness  and  constraint  matrices  are 
employed  in  the  following  arguments,  as  well  as  the  stipulation  that  the  damping  matrix  is  positive 
semi-definite 


(tx,  x)2  0 
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With  the  identical  change  of  coordinates  employed  in  the  undamped  equations,  one  can  consider 
the  system  defined  below  without  loss  of  generality: 

.-1/2.  „  -1/2 
C  =  M  CM 


X  +  CX  +  KX  =  ^  \  +  F 
Ldx  J 


rd<Di 

Lad*  ■  ° 


The  new  system  of  reaction-free  equations  now  contain  a  term  representing  viscous  damping 


X+  (/-  P)CX+  (/  -P)KX  =  (/  -P)F 


X  +  QCX  +  QKX  =  QF 


(3.2)  Damped  Penalty  Equations 

The  penalty  formulation  employed  in  the  case  of  a  damped,  linear  MDOF  system  has  been  selected 
as  in  [Bayo]  and  fKurdilaj.  In  addition  to  the  penalized  kinetic  and  potential  energies,  a  generalized 
Rayleight  dissipation  function  is  introduced 


It-  1  6  B<i> 
T,  =  ix'x+i — — 

£  2  2  e 

1  T  l  <J>TB<J> 
V,  =  ^XTKX+~ — — 
e  2  2  e 
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X£  +  (2(e)CX£  +  2(e)A:XE  +  H/>(E)XE  +  a/,(E)XE  =  (2(e)  F 

(3.3)  Boundedness  of  the  Solutions  of  the  Damped  Penalty  Equations 

Again  before  estimating  the  error  in  the  approximation,  as  in  section  (2),  it  is  necessary  to  establish 
that  the  solution  of  the  penaJty  equations  corresponding  to  damped  system  is  bounded 

FeWISCj 
ll*e(0  | 

for  any  arbitrary,  finite  time  0<t<T. 


HO 


Taking  the  inner  product  of  the  damped  penalty  equations  and  the  derivative  of  the  penalized  so¬ 
lution,  one  can  write  a  time-rate-of-change  of  energy  as  carried  out  in  the  undamped  analysis. 

w  l  t  rd<J?iTrdd>i  .  a  T 

5r{ 2iXe ( Lax]  LwJ)X£  +  2EX£(Lax]  [ax])X£}  + 

T 

Li  .  t  rdOi 

<CX£,Xc>+gx[([^J  ^J)X£  =  <F,X£) 


Because  the  damping  matrix  C  is  symmetric,  positive  semi-definite,  one  can  write  the  inequality 

d  r  1 ;!  v  li  .  tw  v 


^{^!X£|j  +  <XX£,X£>}  + 


^  i  . t  '90^rr3<I>^  a  T  r0dnrr9<Dn 

:5xJ)Xs+2iX‘(L5xj  J*]**.*  *<**«> 


Following  precisely  the  same  steps  as  in  the  analysis  of  the  undamped  case,  one  integration  yields 


1|X£(0  ||2+i<XX£(r),X£(r)> 

*  {^X£(0)||2+i(A:X£(0),X£(0)>+lj'|!F(x)||2rfx} 
+  U‘Q  (|iX£(x)||2+(XX£(x),X£(x)))dx 


\ !!*c  (0  'i2-^  \  (XX  t  (/) , X£  (t) )  £  g  (t)  +  i  J'0  (||X£  (x)  f  +  i  {KXt  (x) ,  X£  (x) »  dx 
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which  enables  one  to  apply  Gronwall’s  inequality  to  conclude 


||*e(0  1  ^  1C, 

<tfX£(r),X£(f)><  k2 


for  a  fixed  time  0<t<T. 

(3.4)  Variational  Dependence  of  Error  on  Constraint  Violation 

The  final  variational  relationship  between  the  constraint  violation  and  the  approximation  error  can 
now  be  achieved  by  subtracting  the  exact,  reactionless  equations  and  the  penalized  equations. 

X-Xt  +  QCX-Q(e)CXt 
+  QKX-Q(e)KXf. 

=  «2-e(e))F  +  pP(e)XE  +  aP(e)XE 

As  in  the  analysis  of  the  undamped  equations,  the  addition  and  subtraction  of  identical  terms  and 
the  application  of  the  triangle  inequality  enables  one  to  write 

*  -  *  til  +  li  Q  -  Q  (e)  II 1  CXE||  +  ii  (211 1|  C  (X  -  XE)  1 
+  !i(2-(2(e)||||XX£||4||Cll||X(X-X£)|| 

<  Q-Q{e)  IIIIFII  +  a||PXE|]  +  a||P (e)  - P||||X£|| 

+  n:jPXEj|+n!|P(e)-P||||XE|| 


In  the  limit  as  e  approaches  0,  one  can  use  the  fact  that 

G(e)  ~>Q 
P(e)  ->P 
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in  operator  norm  and  the  boundedness  of  Xe  to  obtain 


(4.0)  CONSERVATIVE  FORCING  FUNCTIONS 

The  variational  error  bounds  discussed  so  far  have  special  significance  when  the  forcing  funtion  F 
is  in  fact  conservative.  When  this  is  the  case,  one  can  write  [Kurdila] 


Ee(t)  =  E( 0)  -J{(CXE,XE>+^x[ 
o 


-ao>- 

1  r9<b' 

_ax_ 

Lax. 

)X£}dt 


where 


E£  (t)  =  ^xImx£+^xt£kxe 

+  J_  <hrd>  +J-C>racI> 

2e  e  £  2e  E  £ 

E  (0)  =  1xJa/x0  +  ^xJxx0 


As  shown  in  [Kurdial  1,2,31,  one  has  the  bounds 


[<i)E|  +  [|<DE|j  <  kxeE  (0) 


for  some  constant  ic.  This  inequality  implies  that 


l 

;j  c  (X  - xt)  D  +  j| nr  (X  - xt)  |  <K2(e£(0))2 


which  provides  a  rate  of  convergence  when  the  forcing  term  is  conservative. 
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(5.0)  CONCLUSIONS 


By  using  standard  analyses  from  the  field  of  partial  differential  equations  and  linear  regression,  in¬ 
novative  variational  estimates  have  been  derived  that  allow  one  to  monitor  the  fidelity  of  penalty 
method  approximations.  The  analysis  herein  is  applied  to  linear  systems  and  guarantees  conver¬ 
gence  rates  for  conservative  forcing  functions.  The  results  provide  a  qualitative  distinction  between 
the  linear  and  nonlinear  cases:  convergence  of  constraint  norm  to  zero  implies  convergence  of  the 
method  in  the  linear  case,  while  bifurcation  phenomenon  preclude  a  similar  conclusion  in  nonlinear 
simulations,  unless  additional  assumptions  regarding  regularity  are  made. 
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